Abstract
Cutaneous leishmaniasis (CL) caused by Leishmania braziliensis results in chronic skin ulceration and remains challenging to treat. While human transcriptomic studies have identified pathways driving immunopathology, the early events of infection and the molecular transitions from lesion formation to healing are still poorly understood. Here, we performed a longitudinal transcriptomic analysis of skin lesions and draining lymph nodes (dLNs) in the BALB/c ear dermal model infected with L. braziliensis, which recapitulates features of human CL. Using bulk RNA sequencing at 2, 6, and 48 hours and at 14, 35, and 77 days post-infection, we characterized differential gene expression, pathway enrichment, and gene co-expression networks. Ulcerated mouse lesions (Day 35) recapitulated 77% of the inflammatory pathways described in human CL, with many persisting at Day 77 despite “clinical healing”. Mice displayed additional upregulation of genes linked to macrophage polarization (Il12a, Il12b, Il4), nitric oxide metabolism (Arg1, Nos2) and epidermal differentiation (e.g., Crnn, Rptn, Tchh, Lce members). Gene co-expression analysis revealed stage-specific gene modules (M) associated with early innate responses (M3), tissue damage (M1), epithelial-mesenchymal transition (M4), and skin barrier remodeling (M6). A long non-coding RNA-enriched module (M2) was selectively downregulated during the ulceration. Cross-species comparison of ulcerated lesions revealed 16 conserved microRNAs and 12 shared epigenetic regulators, including Mir155, Mir142, Sp140, and Kdm6b, with known roles in inflammation and tissue repair, representing promising host-directed therapeutic targets. Together, this study provides a comprehensive temporal framework of host responses to L. braziliensis and identifies actionable non-coding RNAs and epigenetic pathways with translational potential for CL therapy.
1 Introduction
Leishmaniasis comprises a spectrum of diseases caused by more than 20 species of intracellular protozoan parasites of the genus Leishmania. Clinical manifestations vary widely depending on parasite species and host immune status (1). Cutaneous leishmaniasis (CL) is the most common form and treatment efficacy remains suboptimal, with cure rates as low as 50% (2). In the Americas, Leishmania braziliensis is the principal etiological agent and a major cause of chronic, ulcerative disease (3, 4).
Advances in omics technologies have enabled systems-level investigation of host-Leishmania interactions. Transcriptomic studies of human L. braziliensis lesions have revealed central immunopathological pathways. Novais et al. described a hypothetical metapathway leading from CD8+ T cell activation and cytolysis to IL-1β production, which drives immune-mediated tissue damage (5). To understand the progression of L. braziliensis disease, Christensen et al. used RNA sequencing (RNA-seq) to profile early and late human lesions, demonstrating that these clinical stages were not transcriptionally distinct (6). Additional work showed that chronic, treatment-refractory lesions exhibit increased activation of components of the cytolytic machinery (GZMB, PRF1, GNLY) (7). Across these studies, strong interferon-stimulated gene (ISG) and cytolytic programs consistently emerged (5–7) and were later detected systemically in peripheral blood of L. braziliensis-infected individuals (8). More recently, alterations in the skin microbiota have also been linked to disease severity (9).
Murine models have provided complementary mechanistic insights. Single-cell RNA sequencing (scRNA-seq) of L. major infection in C57BL/6 mice demonstrated macrophage remodeling within ear lesions, along with downregulation of initiation factor-2 (EIF2) and several ribosomal subunits (10, 11). Another scRNA-seq study demonstrated cellular remodeling caused by L. major, suggesting that dermis-resident macrophages orchestrate Th2-type immune responses along with innate lymphoid cells and eosinophils within the lesion, serving as a replicative niche for L. major during the progression of more severe CL (12).
Despite these advances, our understanding of the temporal sequence of events during in vivo infection remains limited, particularly with respect to early host responses, the transition to ulceration, and the dynamics of lesion healing. Furthermore, gaps remain in understanding the communication between the immune responses occurring in the skin and the draining lymph node (dLN). Although the BALB/c ear dermal model of L. braziliensis infection has been widely used for nearly two decades—for example, in preclinical vaccine and drug screening studies—and recapitulates key features of human CL, including ulceration, parasite dissemination to lymphoid tissues, and a dominant Th1-type response (13), no longitudinal, genome-wide transcriptomic analysis has been performed in this model.
In parallel, there is growing interest in host-directed therapies (HDTs) to improve clinical outcomes in leishmaniasis. In cutaneous leishmaniasis, particularly in Leishmania braziliensis infection, tissue damage and treatment failure are often driven by dysregulated host inflammatory responses rather than parasite burden alone. HDTs aim to modulate host immune, metabolic, or epigenetic pathways rather than directly targeting the parasite, with the potential to mitigate immunopathology, reduce lesion chronicity, and enhance the efficacy of antiparasitic drugs (14, 15). Epigenetic regulators, microRNAs, and inflammatory signaling pathways have emerged as promising intervention points, yet few have been evaluated in vivo in the context of L. braziliensis infection.
In this study, we conducted a comprehensive time-course RNA-seq analysis of L. braziliensis infection in the murine ear, profiling both lesions and matched dLNs from early infection through ulceration and healing. By integrating differential gene expression, pathway enrichment, co-expression networks, and cell-type deconvolution, and by comparing our data with multiple human and murine CL datasets (5–7, 9–11), we define conserved and species-specific molecular programs underlying disease progression. We also identify non-coding RNAs and epigenetic regulators with translational potential as host-directed therapeutic targets.
2 Materials and methods
2.1 Ethics statement
All procedures involving animals were approved by the Ethics Committee on the Use of Animals (Comissão de Ética no Uso de Animais – CEUA) of the Instituto Gonçalo Moniz (IGM-FIOCRUZ, Salvador, Bahia-Brazil) under license no. 008/2014 and complied with national regulations for laboratory animal welfare. All procedures involving euthanasia were performed in accordance with institutional and national guidelines for the humane treatment of laboratory animals.
2.2 Parasites and culture
L. (Viannia) braziliensis (MHOM/BR/01/BA788) virulence was maintained through passages in BALB/c mice. Parasites were grown in vitro until reaching stationary phase in a B.O.D. incubator at 26 °C in Schneider’s insect medium (Invitrogen) supplemented with 10% inactive Fetal Bovine Serum (FBS) (Gibco), 2 mM L-glutamine, 100 U/ml penicillin and 100 µg/ml streptomycin (all from Invitrogen) (13). Stationary-phase promastigotes were used for infection.
2.3 Animals and experimental design
Twenty-one female BALB/c mice (6–8 weeks old), were obtained from IGM-FIOCRUZ animal facility. A total of 18 mice were inoculated intradermally in the center of ear pinna with 105 stationary-phase L. braziliensis promastigotes in 10 μL of saline using a 27.5-g needle (13); three non-manipulated mice were kept as controls. Ear thickness was recorded weekly with a digital caliper (Thomas Scientific). Mice (n=3 per time point) were euthanized at 2 h, 6 h, 48 h, 14 days, 35 days and 77 days post-infection (dpi) (Figure 1). Animals were euthanized using carbon dioxide (CO2) inhalation in a dedicated euthanasia chamber with a gradual fill rate of 20% of the chamber volume per minute, with the CO2 flow maintained for at least one minute after clinical death, followed by confirmation of death prior to tissue collection. For each animal, the infected ear was collected and immediately snapped frozen and stored in liquid nitrogen until RNA extraction. The retromaxillar lymph nodes draining this ear (dLNs) were collected, fixed overnight in RNAlater solution (Ambion) and then stored in liquid nitrogen until RNA extraction. For clarity, lesions were classified into an early phase (2 h–48 h post-infection, without overt clinical signs), a pre-ulcerative phase (Day 14, characterized by papule formation without tissue breakdown), and a late phase comprising the ulcerated stage (Day 35) and the healed stage (Day 77).
Figure 1
2.4 RNA extraction
The cryopreserved ears were processed individually into a fine powder using a porcelain mortar/pestle in the presence of liquid nitrogen. Powdered tissue was homogenized in RLT buffer (Qiagen) for total RNA isolation using the RNeasy Mini kit (Qiagen) following manufacturer’s instructions. To avoid RNA degradation at later time points, tissue was thawed only in the presence of RLT and vigorously vortexed. dLNs were extracted with the RNeasy Micro Kit (Qiagen) following the manufacturer’s instructions. RNA quantity and purity were assessed by NanoDrop 1000 (Thermo Scientific, Waltham, MA) and RNA integrity by Agilent Bioanalyzer (Agilent Technologies, Santa Clara, CA). RNA was frozen at −80 °C until library construction.
2.5 Library preparation and sequencing
For ear samples, ribosomal RNA was depleted using the Eukaryote RiboMinus kit (Invitrogen,Carlsbad, CA). Strand-specific libraries were prepared using TruSeq Stranded Total RNA Library Prep workflow (Illumina, San Diego, CA) according to the manufacturer’s recommendations. Libraries were sequenced on an Illumina NextSeq500 system using v2 chemistry in High Output mode, configured for 150 sequencing cycles, generating 75 base pair paired-end reads. For dLNs samples, libraries were prepared using the TruSeq Stranded mRNA Library Prep Kit (Illumina, San Diego, CA) and sequenced on an Illumina HiSeq 2500 system using a stranded paired-end protocol, generating 75 base pair paired-end reads. Target sequencing depth were ~30 M reads per ear sample and ~20 M reads per dLN sample. Library QC metrics and per-sample read counts are provided in Supplementary Table 1.
2.6 Gene expression data analysis
Read quality control (QC) was assessed using FastQC (http://www.bioinformatics.babraham.ac.uk/projects/fastqc/). Low-quality bases and adapter sequences were removed with Trimmomatic v0.36 (21). Reads passing QC were aligned to the mouse reference genome from GENCODE release M21 (assembly GRCm38.p6) (22) using STAR v2.7.1a with default parameters (23). Host-unmapped reads were subsequently re-mapped to the L. braziliensis M2904 genome (24). Annotation (GFF) and sequence (FASTA) files were obtained from the Leish-ESP database (http://leish-esp.cbm.uam.es/). Transcripts counts were generated using featureCounts. The Trimmed Mean of M-values (TMM) normalization was performed to obtain normalization factors across samples and to account for underlying RNA composition (25). Prior to testing for differential expression, low-expression genes were filtered using the default parameters of the filterByExpr function implemented in edgeR v. 3.26.8 (26). P-values were corrected for multiple testing using the False Discovery Rate (FDR). Genes were considered differentially expressed (DEGs) when FDR ≤ 0.05 and absolute log2 fold change ≥ 2 for lesion samples or ≥ 1.5 for dLN samples, unless otherwise specified. For comparative analyses, we also incorporated publicly available gene expression datasets from human cutaneous leishmaniasis (CL) lesions, referred to as Hs (5), Hs2 (6), Hs3 (7), and Hs4 (9).
2.7 Gene set enrichment analysis
To assess genome-wide expression patterns and identify coordinated changes in predefined gene sets across different time points, we performed Gene Set Enrichment Analysis (GSEA, Broad Institute; https://www.broadinstitute.org/) (27). For each time point, genes were ranked by their differential expression p-values, and enrichment was evaluated against curated gene sets. GSEA procedures were carried out as previously described (28).
In addition to our dataset, we applied GSEA to publicly available expression profiles from human CL lesions (Hs) (5), L. major-induced mouse lesions (Lm) (10), and popliteal lymph node aspirates from L. infantum–infected dogs (Li) (16). To account for nonspecific responses resulting from needle injury and inoculation, we also included a dataset from saline-injected mouse ears (needle control, NC) (17). Gene sets used in the analysis consisted of curated canonical pathways from the Kyoto Encyclopedia of Genes and Genomes (KEGG) (29) and the Reactome Pathways Database (30).
2.8 Estimation of cell composition
Digital cytometry was performed to characterize the immune infiltrate using the seq-ImmuCC algorithm, which estimates 10 immune cell populations from RNA-seq count data. The method applies least-squares linear regression machine-learning approach to reference signature matrices specifically developed for different mouse tissues. For our analyses, we used the PBMC matrix for skin lesions and the spleen matrix for dLNs (31). Visualizations were generated in R version 3.4.4 (http://www.r-project.org/) using base functions and the ggplot2 package (32). Statistics between groups were assessed by Mann-Whitney U-test using the rstatix package.
To estimate the cell composition from the gene modules identified by CEMiTool (see below), we used the Web-based Cell-type Specific Enrichment Analysis of Genes (WebCSEA) (33) (https://bioinfo.uth.edu/webcsea/). This resource encompasses 111 scRNA-seq panels of human tissues and 1,355 tissue-cell types from 61 tissues across 11 human organ systems. Mouse module genes were first converted to their human orthologs using g:Profiler (https://biit.cs.ut.ee/gprofiler/orth). Enrichment analysis was then performed on WebCSEA, focusing on the top 20 general cell-type datasets. Only cell types with enrichment above the Bonferroni-corrected significance threshold of 3.69 × 10-5 (p-value set at α = 0.05, adjusted for 1,355 comparisons) were displayed.
To identify macrophage subtypes and other cellular populations in these samples—and to compare them with subtypes predicted from bulk RNA-seq data of human CL biopsies (9) - we applied deconvolution techniques using single-cell RNA-seq datasets from murine and human models of acute skin wounds induced by punch biopsy (18–20). A detailed description of all methods can be found in Supplementary Methods, Appendix.
2.9 Co-expression network and enrichment analysis
To identify groups (i.e., modules) of co-expressed genes in our dataset, we used the Co-Expression Modules identification Tools (CEMiTool) R package, which enables automated module detection through an unsupervised gene filtering method coupled with automated parameter selection and enrichment analysis. Data used in this analysis was normalized as previously described (34). CEMiTool also performs over-representation analysis (ORA) to identify enriched pathways and constructs interaction networks to detect protein-protein interaction (PPI) patterns.
For pathway enrichment, we used the Hallmark gene sets from Mouse Molecular Signatures Database (MSigDB). To identify protein interactions relevant to innate immune responses, we queried the InnateDB database (35). Fisher’s exact test was used to calculate the overlap p-value for each analytic tool; significance was attributed to p‐values < 0.05.
3 Results
3.1 Gene expression dynamics reveal temporal phases of lesion development and resolution
In the dermal infection model with L. braziliensis, ear thickness begins to increase around Day 14 after inoculation, with lesion volume peaking between Days 21 and 35, followed by spontaneous healing around Day 70 (Figure 1). Accordingly, sampling time points were selected to capture key phases of infection and lesion evolution, including early host responses, disease progression, ulceration, and resolution. Specifically, early time points (2 h and 6 h) were chosen to reflect immediate innate responses and initial parasite–host interactions, while 48 h captures the transition toward early adaptive immunity phase preceding macroscopic lesion development. Sampling at Day 14 corresponds to the late pre-ulcerative stage, characterized by papule formation. Day 35 represents peak lesion size with ulcerated morphology, whereas Day 77 corresponds to the post-healing phase, marked by lesion resolution. For analytical purposes, we defined the early phase as 2 hours to 48 hours, Day 14 as pre-ulcerative phase, and the late phase as Day 35 (ulcerated lesions) and Day 77 (healed lesions).
A total of 20 ear samples and 21 dLN samples were collected across the following groups: Control, 2h, 6h, 48h, Day 14, Day 35, and Day 77, and analyzed by RNA-seq. Principal component analysis showed that ear samples clustered tightly by time point, indicating a clear temporal pattern in gene expression. Notably, the 48 h time point clustered more closely with control samples (Supplementary Figure 1A, Appendix). In contrast, dLN samples exhibited greater variability among biological replicates. Even so, Day 35, corresponding to peak lesion development, showed a distinct transcriptional profile compared to other time points (Supplementary Figure 1B, Appendix).
While the sequencing depth was not sufficient for quantitative analysis of Leishmania transcript abundance, the proportion of reads mapping to the L. braziliensis genome suggests that parasite transcriptional activity was highest at Day 35 in both ear lesions and dLNs (Supplementary Table 1, Supplementary Figures 1C, 1D, Appendix).
3.2 Lesion analysis
3.2.1 The dermal model recapitulates key inflammatory gene expression patterns observed in human CL biopsies
To investigate transcriptomic signatures associated with the course of infection, we performed differential gene expression (DEG) analyses comparing each infected time point to the control group. The transcriptomic profile revealed a peak in DEGs at Day 35 (1,817; 1,017 up, 800 down), followed by Day 77 (1,111; 990 up, 121 down), 6 h (390; 366 up, 24 down), Day 14 (212; 208 up, 14 down), 2 h (207; 183 up, 24 down), and 48 h (70; 69 up, 1 down), with most genes being upregulated (Figure 2A). Volcano plots illustrate the magnitude of gene expression changes in infected ears (Supplementary Figure 2, Appendix). All DEGs are listed in Supplementary Table 2.
Figure 2
Noteworthy, several DEGs identified in the murine model correspond to genes proposed by Novais and colleagues as part of a “metapathway” driving immunopathology in human CL (5), highlighted in blue in Figure 2B and in Supplementary Figure 2, Appendix. These include immune proteasome components (Ifng, Tnfa and Psmb9),chemokines involved in cell recruitment (Cxcl10, Cxcl9, Cxcl5, Cxcl2, and Cxcl1),cytotoxicity-associated genes (Gzmb, Prf1 and Casp4), and inflammasome-related genes (Nlrp3, Aim2, Casp1, Il1b, and Gbp5). Additional genes from this “metapathway” (Psma4, Psme2, Casp3, Casp7, and Bid) were also regulated but below the fold-change threshold to be classified as DEGs (Supplementary Table 3).
Unexpectedly, most genes within this “metapathway” remained upregulated even in thehealed lesions at Day 77. Furthermore, we did not observe significant inhibition of the seven genesassociated with the skin barrier pathway previously described by Novais et al. (5) (Supplementary Table 3).
3.2.2 Divergent immune and skin barrier gene expression between human and mouse lesions
To further investigate the expression of skin barrier-related genes, we expanded our analysis to include key genes identified by Christensen et al., 2016 (Hs2), who reported the regulation of antimicrobial and skin barrier-associated genes in human CL biopsies. When comparing transcriptional profiles from mice at Days 35 and 77 with those from humans (Hs2), we observed a notable overlap of B cell markers (Cd79a, Cd19, and Cd20 [Ms4a1]) between Day D77 and Hs2 (Figure 2C). Inflammatory, inhibitory, and activation markers showed broadly similar expression profiles between mice and humans.
In contrast, mouse lesions at Days 35 and 77 exhibited marked upregulation of key genes involved in macrophage polarization (Il12a, Il12b, and Il4), as well as genes associated with L-arginine metabolism and nitric oxide production (Arg1, Nos2 and Cat-2 [Slc7a2]) (Figure 2C). Moreover, unlike the human biopsies, where these genes were significantly downregulated relative to healthy skin, mouse lesions showed no evidence of downregulation in skin barrier genes (e.g. Flg, Flg2, Lor, Dsc1) when compared to non-manipulated controls. Instead, multiple genes from the epidermal differentiation complex (EDC), including Crnn, Tchh, Tchhl1, Rptn, Lce3a, Lce3d, Lce3e, Lce3c and Lce3b, were significantly upregulated in mouse lesions (Figure 2C). A quantitative summary highlighting the contrasting directionality of EDC regulation between species is provided in Supplementary Figure 3, Appendix. These findings are further supported by complementary analyses of Krt and Krtp gene regulation, which are presented in Section 3.2.7 (Gene co-expression analysis during lesion progression and resolution).
3.2.3 The ulcerated and healed phases exhibit the greatest overlap in differentially expressed genes
To identify transcriptional responses associated with each phase of lesion development and the transitions between them, we generated UpSet diagrams from the list of DEGs (both up- and down-regulated) across the sampled time points. This analysis revealed that the largest overlap occurred between Day 35 and Day 77, comprising 406 shared upregulated genes and 57 shared downregulated genes (Supplementary Figures 4A, B, Appendix).
We then examined the top 20 DEGs shared across at least two experimental time points. Among the upregulated genes in the lesions, some immunoglobulins transcripts were highly expressed at both Day 35 and Day 77 (Supplementary Figure 4C, Appendix). Moreover, key chemokines involved in cell recruitment, such as Cxcl5, Cxcl3, Cxcl2, and Ccl4, were identified in the intersection of 2h, 6h, Day 35, and Day 77. Notably, several genes previously implicated in leishmaniasis, including Fpr1, Trem1, Acod1, Gzma, and Olfm4, were also prominently represented. We further identified four strongly upregulated predicted genes: Gm4841 (an interferon-inducible GTPase), Gm15056 (a defensin), and Gm11555 and Gm7544 (both keratin-associated proteins).
In contrast, the top 20 downregulated transcripts showed far less overlap across time points, and unlike the upregulated genes, they could not be readily grouped into functional categories. Interestingly, 13 of the downregulated DEGs were computationally predicted genes annotated with the “Gm”prefix, highlighting a predominance of poorly characterized genes among the downregulated DEGs (Supplementary Figure 4D, Appendix).
3.2.4 GSEA reveals distinct temporal phases and species-specific immune and epidermal signatures in L. braziliensis lesions
To gain insight into gene regulation, even at time points with few DEGs, we initially performed KEGG pathway enrichment analyses using the full set of genes detected by RNA-seq in lesion samples. To account for nonspecific responses caused by mechanical stress from needle injury and inoculation, we included a control dataset from mouse ears injected with saline (needle control, NC) (17). We also incorporated datasets from human lesions caused by L. braziliensis (5) (Hs) and from L. major-induced CL in C57/BL6 mice, early sampled 28 days post-infection (10) (Lm). Following manual curation and removal of redundant or overlapping gene sets, a heatmap based on normalized enrichment scores (NES) revealed gene sets expressed exclusively at the early phase of lesion development, those shared between early and late phases, and others restricted to the late phase (Figure 3A). The pre-ulcerative phase at Day 14 displayed a distinct profile.
Figure 3
In the exclusive early response (labeled in grey), ECM receptor, VEGF, and MAPK signaling pathways were activated as early as 2 hours post-infection (Figure 3A). Several additional pathways triggered during the early phase (2-48h) were also reactivated later, at Days 35 and 77 (labeled in green), but with greater intensity. These included pathways involved in host defense, immune cell activation, and cell recruitment (e.g., TLR signaling, NK cell cytotoxicity, apoptosis, Leishmania infection, leukocyte migration, and B and T cell signaling).
Late-phase pathways (Days 35 and 77) that overlapped with human lesion biopsies (labeled inpurple) were enriched for T cell cytotoxicity, mast-cell activation, suppression of amino acid andfatty acid metabolism, and decreased melanogenesis. Notably, some host defense-related pathways appeared only in late-stage mouse lesions and were not detected in human samples. These included epithelial signaling in Helicobacter pylori infection, Escherichia coli pathogenicity, proteasome activity, and arginine and proline metabolism. All KEGG-identified pathways are listed in Supplementary Table 4 and the number of pathways detected at each time point is shown in Figure 3B.
To complement the KEGG analysis and reveal additional functional differences, we performed Reactome pathway enrichment analysis. The most striking observation was a broad transcriptional suppression during the pre-ulcerative phase at Day 14 (Figure 3C, Supplementary Table 5). Although lesions were papular and non-ulcerated at this time point, 88 Reactome and KEGG pathways were downregulated, spanning extracellular matrix remodeling, immune and inflammatory responses, cell signaling and migration, as well as metabolic and mitochondrial programs (Figures 3B, C). In total, 92 pathways were grouped into seven central biological programs, as detailed inSupplementary Table 5 (Classification of pathways into central biological programs).
Notably, ribosomal and translational pathways were also downregulated, suggesting that reduced biosynthetic capacity may underlie this pleiotropic transcriptional repression. In contrast, only four pathways were upregulated, all related to epidermal differentiation and skin barrier maturation. Accordingly, we refer to Day 14 as a pre-ulcerative (transcriptionally silent) phase, reflecting widespread pathway-level suppression despite the presence of early clinical manifestations.
We next examined pathways showing opposite enrichment patterns between mouse ulcerated lesions(Day 35) and human CL biopsies. Among the 351 pathways analyzed, only two displayed contrastingregulation: Formation of the Cornified Envelope and Keratinization, both upregulated in mouse lesions but downregulated in human samples. Notably, both pathways were among the four already upregulated during the pre-ulcerative (transcriptionally silent) phase at Day 14, highlighting early activation of epidermal differentiation programs in the murine model. All identified pathways are listed in Supplementary Table 5.
This analysis highlights four key findings: 1. A distinct pre-ulcerative “transcriptionally silent” phase at Day 14, characterized by marked downregulation of immune, metabolic, signaling, and protein synthesis-related pathways (e.g., Ribosome, Complement, and Focal Adhesion in KEGG), alongside selective upregulation of epidermal development and skin barrier pathways (Reactome); 2. L. major-induced lesions at Day 28 (Lm) represent a transitional state, with enrichment of NK-cell cytotoxicity, Complement, and Leishmania infection pathways; 3. Pathway profiles at Days 35 and 77 were highly similar; and 4. Although mouse and human L. braziliensis lesions shared many enriched pathways during the ulcerated phase, mouse lesions showed additional activation of host-defense pathways and consistent upregulation of epidermal development and barrier function pathways—processes that were instead downregulated in human CL.
3.2.5 Macrophages and natural killer cells predominate during the ulcerated phase of the lesion development
To estimate immune cell-type composition across lesion time points, we applied CIBERSORTx deconvolution to the bulk RNA-seq data using the ImmuCC PBMC signature matrix, given the absence of a skin-specific matrix. Our initial rationale for using PBMC-derived reference matrix to analyze skin transcriptomic data was to capture circulating immune populations recruited to the skin following infection, while a spleen-derived matrix was the closest available proxy for secondary lymphoid tissue. Most samples showed a predominance of monocytes, likely reflecting the use of a PBMC-derived reference. Significant shifts in immune composition were observed only at Days 35 and 77 (Figure 4A).
Figure 4
At Day 35, macrophage abundance increased markedly while monocytes decreased. The sharp rise in natural killer (NK) cells was particularly notable (Figures 4A, B). By Day 77, macrophages remained elevated, and dendritic cells also showed a substantial increase (Figures 4A, B).
3.2.6 Murine and human lesions exhibit convergent macrophage dynamics and distinct epidermal responses
To further characterize macrophage subtypes and compare cell populations with those in human CL biopsies (9), we performed computational reference-based deconvolution using single-cell RNA-seq datasets derived from murine (Figure 4C) and human (Figure 4E) acute cutaneous wounds (ACW) generated by punch biopsy. At the time these analyses were performed, no publicly available single-cell RNA-seq datasets from human cutaneous leishmaniasis lesions were available. We therefore opted to use reference datasets derived from acute skin inflammation as a pragmatic approach to enable cross-species comparative analyses. Heatmaps show the estimated proportions of major cell clusters in each sample, inferred using the MuSiC algorithm and the corresponding ACW reference datasets (Figures 4D, F). As these reference datasets are derived from acute wound models, the inferred cell-type proportions likely reflect relative functional states rather than absolute cellular identities in chronic leishmaniasis lesions.
In murine lesions, M2-like (regulatory/tissue-adaptive) macrophages were more abundant at Day 35, while M1-like (inflammatory/antigen-processing) macrophages predominated at Day 77 (Figures 4D, G). Human CL lesions similarly exhibited a dominant M2-like macrophage profile (Figures 4F, G). The M1-like/M2-like ratio for both species is shown in Figure 4H; however, the estimated balance between macrophage states may be influenced by differences between acute wound reference datasets and chronic infectious lesions.
Keratinocyte abundance remained relatively stable in mouse lesions compared to controls, with a slight upward trend. In contrast, human CL lesions displayed a clear and substantial reduction in keratinocyte proportions (Figure 4I). The full list of identified cell types, estimated cellular proportions, and definingmarker genes is provided in Supplementary Table 6.
3.2.7 Gene co-expression analysis during lesion progression and resolution
GSEA evaluates predefined gene sets based on prior biological knowledge (36). In contrast, gene co-expression network analysis identifies modules of genes that share similar expression patterns, allowing phenotype associations without relying on prior assumptions (37). Using CEMiTool, we identified co-expression modules that characterize different stages of lesion progression and resolution.
Among the eight modules identified, six showed marked variation across the phases of lesion development (Figure 5, Supplementary Table 7). Over-representation analysis (ORA) revealed that two modules were associated with tissue remodeling (M4 and M6), three with immune pathways (M1, M3 and M5) and one module (M2) could not be assigned any significant process. An interactive exploration of the modules, including interaction networks, is available online (38). To further investigate their functional relevance, we used WebCSEA to infer the predominant cell types associated with each module.
Figure 5
Module M6, associated with apical junction pathways, was the earliest module to show activation. It peaked at 6 hours post-infection and again at Days 14 and 35, coinciding with papule formation and ulcer opening, respectively (Figure 5). M6 was enriched for keratin-related genes: 37% of its members were keratins (KRTs) or keratin-associated proteins (KRTAPs). Of 79 co-expressed Krt/Krtap genes, 74 were grouped in M6 and were all upregulated during infection (Supplementary Figure 5A, Appendix). In contrast, human lesions showed broad downregulation of these genes, with 81% of 43 regulated KRT/KRTAP genes being repressed (Supplementary Figure 5B, Appendix). Given the role of keratins in skin structure and repair, this module likely reflects epithelial remodeling during early infection, ulcer formation, and healing.
Module M4, also active during the early stages preceding ulcer opening, was enriched for epithelial–mesenchymal transition and apical junction genes. It includes collagens and adhesion-related genes whose expression declined progressively from early infection through Days 14 and 35 (Figure 5). WebCSEA linked M4 to mesenchymal lineages, including fibroblasts, chondrocytes, stromal cells, and other mesenchymal populations, consistent with its structural and matrix-related functions (Supplementary Figures 6A, B, Appendix).
Among immune-related modules, M1 and M3 showed similar patterns, with strong upregulation during the ulcerated phase (Days 35–77). M3 also displayed moderate activation early in infection (2–6 h).
Module M1 showed a classic IFN-γ signature, including genes such as Irf1, Irf7, Irf8, Dhx58, Isg15, and Stat1, along with genes linked to antigen presentation (Ctss), tissue damage (Gzmb), and inflammasome activity (Casp1) (Supplementary Figure 6C, Appendix). WebCSEA linked M1 to lymphoid populations, including T cells, NK cells, and dendritic cells, highlighting its association with adaptive immune response (Supplementary Figure 6D, Appendix).
M3 interaction network, in turn, included key regulators of immune cell recruitment and activation (Cd14, Cd40, Vav1, Itgb2, Ccr1, Ccr5, Cxcr2, Ccl2, Nlrp3) along with cytokine signaling genes (Socs3, Fgr), reflecting its role in innate immunity (Supplementary Figure 6E, Appendix). WebCSEA linked M3 to monocytes, macrophages, neutrophils, dendritic cells, and myeloid cells (Supplementary Figure 6F, Appendix).
Module M5 displayed a distinct late-phase pattern, with marked upregulation at Day 77. WebCSEA linked this module to B cells (Supplementary Figure 7A, Appendix). M5 contained 66% immunoglobulin (Ig) genes. Early isotypes (Ighm, Ighd, and Ighg1) were already expressed at Day 35, whereas additional isotypes (Igha, Ighe, Ighg2b, Ighg2c, Ighg3) were upregulated at Day 77 (Supplementary Figure 7B, Appendix). A heatmap of Ig-related DEGs confirmed strong expression in lesion tissues (Supplementary Figure 7D, Appendix), and peak expression in dLNs at Day 35 for several isotypes (Supplementary Figure 7C, Appendix).
Module M2 exhibited a unique pattern, with downregulation at Day 35 followed by recovery at Day 77. It was not enriched for any known pathways or cell-type signatures, and no coherent functional theme emerged from the interaction network analysis (38). Manual annotation using MGI identifiers revealed that M2 consisted primarily of unclassified genes (52%), non-coding RNAs (23%), pseudogenes (10%), and only 16% protein-coding genes (Supplementary Figure 8, Appendix, Supplementary Table 7). Notably, most downregulated genes were ncRNAs, pseudogenes, or unclassified elements, whereas upregulated genes were primarily protein-coding. Among them, Mcpt1, Mcpt2, Mcpt8 (mast-cell proteases) and Arg1 were notable, though most other upregulated genes lacked clear functional annotation and appeared associated with macrophage or cytotoxic T cell activity.
3.3 Draining lymph nodes analysis
To investigate transcriptomic signatures associated with infection dynamics in lymph nodes (dLNs), we applied the same analytical workflow used for lesion samples.
3.3.1 Differentially expressed genes in dLN peak at Day 35 and return to homeostasis by day 77
Gene expression analysis revealed a peak in DEGs on Day 35, with 1,686 genes identified. At this time point, several key components of the immunopathology-associated “metapathway” were highly expressed, including Il1b, Gzmb, Ifng, Cxcl9, Gbp5, Nlrp3, and Cxcl10. In contrast, only modest transcriptional changes were observed at the other time points, with DEGs ranging from 6 at 6 hours to 67 on Day 14. Remarkably, no DEGs were detected on Day 77, suggesting a return to homeostasis (Supplementary Figure 9, Appendix, Supplementary Table 2).
Gene expression overlap analysis showed limited sharing of DEGs across time points. The intersections of [Day 14 with Day 35] and [48 hours with Day 35] exhibited the highest numbers of shared upregulated genes (36 and 35 DEGs, respectively) (Supplementary Figure 10A, Appendix). The [Day 14-Day 35] intersection also showed the largest overlap of downregulated DEGs (Supplementary Figure 10B, Appendix).
We next examined the top 20 DEGs shared across at least two time points. Among these, 12 immunoglobulins genes were upregulated, along with key regulators of B-cell activation and germinal center responses, including Aicda, Nuggc, and Rgs13 (Supplementary Figure 10C, Appendix). Conversely, several genes involved in neutrophils and polymorphonuclear leukocytes responses, such as S100a8, S100a9, Ngp, Mmp8, and Camp, were downregulated, along with genes involved in heme and hemoglobin metabolism (Supplementary Figure 10D, Appendix).
3.3.2 Macrophages and B cells are predominant in the dLN during the ulcerative phase
Immune cell populations in the dLN were inferred using CibersortX, following the same approach used for lesion samples, but employing a spleen-derived signature matrix. The most significant changes occurred on Day 35, characterized by increased proportions of macrophages and B cells, and a concurrent decrease in DC and CD4+ T cells (Supplementary Figure 11, Appendix). These shifts were consistent with the volcano plot expression patterns and with the top 20 shared DEGs (Supplementary Figures 9, 10, Appendix).
3.3.3 Pathway analysis of immune responses in draining lymph nodes
To identify pathways associated with events preceding ulceration and those occurring during healing, we conducted GSEA on the dLN dataset. Unlike lesion tissue, the dLN transcriptome did not exhibit a clear temporal pattern of immune regulation. Nonetheless, GSEA revealed molecular signatures that were not captured by DEGs analysis alone.
Day 35 was characterized by substantial downregulation of pathways related to signal transduction, cellular differentiation, and proliferation. Interestingly, the antigen presentation pathway was activated as early as 2 hours post-infection and remained active until Day 35. In contrast, apoptosis-related pathways were upregulated only during the healing phase (Day 77) (Supplementary Figure 12A, Appendix). Notably, unlike what was observed in lesion tissue, we did not detect a pronounced downregulation of immune-related pathways or protein translation processes in the draining lymph nodes (dLNs) at Day 14.
To contextualize the dLN response, we compared our dataset with a published transcriptomic profile of popliteal lymph nodes aspirate of dogs experimentally infected with Leishmania infantum (Li), a visceralizing species, sampled 28 days post-infection (16). The comparison revealed only six shared pathways (23%): Lysosome, Chemokine signaling, Proteasome, Antigen presentation, NK-mediated cytotoxicity, and Apoptosis (Supplementary Figure 12A, Appendix). Noteworthy, the Leishmania infection pathway was absent from theLi dataset. All KEGG-identified pathways for the dLN are listed in Supplementary Table 4, and the number of pathways detected at each time point is shown in Supplementary Figure 12B, Appendix.
Gene co-expression analysis of the dLN did not reveal temporal signatures or enriched pathways as clearly as in the lesions. An interactive analysis of the modules, including interaction networks, is available online (39).
3.4 Non-coding RNAs and epigenetic enzymes as potential therapeutic targets
We next examined specific classes of molecules recently proposed as potential intervention targets, including non-coding RNAs (ncRNAs) and enzymes involved in epigenetic regulation (Figures 6A, B). The identification of conserved regulated genes and pathways in both mice and humans reinforces the translational relevance of the murine model for the preclinical evaluation of novel interventions strategies.
Figure 6
Although our study did not employ specific miRNA-capture methodologies, we identified 36 miRNAs-associated transcripts and 8 host genes (mir-HGs) in lesion samples, of which 20 were differentially expressed (Figure 6A, Supplementary Table 8). Early regulation was evident as soon as 2 hours post-infection, with most miRNAs/mir-HGs upregulated. A second phase of regulation occurred on Day 35, characterized primarily by downregulation.
Several of these miRNAs have been previously implicated in Leishmania infectionin vitro and are known to play crucial roles in parasite persistence, includingMir221, Mir142, and Mirlet-7e, as well as mature miRNAs derived from Mir205hg, Mir133a-1hg, and Mir155hg (40–45) (Supplementary Table 8). We also identified six miRNAs not previously associated with Leishmania infection: Mir6236, Mir1938, Mir17hg, Mir23a, Mir703, and Mir5119 (Figure 6A). For cross-species comparison, we incorporated human expression data from Amorim et al. (Hs3 and Hs4) (7, 9). Notably, Mir142, Mir100hg, Mir205hg and Mir155hg displayed consistent regulatory patterns across species and datasets, whereas Mir17hg, Mir23b and Mir99a were consistently regulated in at least one dataset (Figure 6A).
In case of epigenetic regulators, we analyzed three enzymes classes (46): 18 DNA-modifying enzymes, 103 histone-modifying enzymes, and 7 proteinarginine methyltransferases (Supplementary Table 9). Among these, 3 DNA-modifying and 9 histone-modifying enzymes were differentially regulated. Of these, 1 DNA-modifying and 8 histone-modifying enzymes displayed consistent regulation in both mouse and human lesions, with most peaking during the ulcerated phase (Figure 6B). Four of these enzymes have available inhibitors or chemical probes, enabling future functional validation.
Tet1, a DNA demethylase that contributes to M1 macrophage polarization, was downregulated both at early (2–6 h) and late (Days 35-77) time points, suggesting it is not a suitable target for boosting anti-inflammatory responses. Similarly, UHRF1, an epigenetic repressor of inflammatory genes in myeloid cells, is unlikely to be an optimal inhibition target. In contrast, enzymes linked to inflammatory processes and cytokine production, such as Sp140 (a bromodomain-containing protein) and Kdm6b (a histone demethylase), were positively regulated, and inhibitors like GSK761 and GSK-J4 warrant further functional evaluation.
4 Discussion
4.1 Mouse lesion closure despite sustained inflammatory response
Several mouse models have been used to dissect the pathogenesis of CL and evaluate potential therapies (47, 48). While lesion kinetics in humans and mice are similar, murine lesions resolve spontaneously within weeks, whereas human healing is markedly slower (13, 49).
To assess how closely the L. braziliensis mouse model reflects human disease, we analyzed temporal changes in DEGs, pathways, and co-expression modules. At Day 35, 83% of genes within a “metapathway” linked to CD8+ T cell activation and IL-1β production were upregulated, mirroring the immunopathology observed in human CL (5). However, key components of the immunoproteasome (Psma4, Psme2) and apoptosis-related genes (Casp3, Casp7, Bid) were less modulated, suggesting less sustained CD8+ T cell activation during the ulcerated phase in mice. The predominance of Gzma over Gzmb expression further diverged from the human profile (50–52). We also found no increase of Prdm1 (Blimp-1) expression in lesions, whereas in C57BL/6 mice, hypoxia in inflamed skin has been shown to induce Blimp-1 during cytolytic CD8+ T cell differentiation (53).
Both the ulcerated (Day 35) and healed (Day 77) phases also showed strong upregulation of genes involved in macrophage polarization (Il12a, Il12b, Il4) and nitric oxide metabolism (Arg1, Nos2) (Figures 2, 7). Despite complete clinical healing by Day 77, many inflammatory pathways remained active, likely reflecting ongoing resolution of inflammation. KEGG analysis revealed a 77% concordance between murine ulcerated lesions and human biopsies. However, mice displayed several additional defense-related pathways and broader macrophage activation signatures, consistent with a more efficient immune response and faster parasite clearance.
Figure 7
Taken together, these data indicate that the mouse model recapitulates key inflammatory features of human CL and mounts a stronger effector response that drives earlier lesion resolution; however, the persistence of inflammatory pathways at Day 77, despite apparent clinical healing, highlights the importance of incorporating molecular endpoints when evaluating therapies in this model.
4.2 Divergent epidermal and wound-healing programs in murine and human cutaneous leishmaniasis
Wound healing is essential for maintaining skin integrity and occurs in three interconnected stages: inflammation, proliferation, and tissue remodeling (54). In L. braziliensis infection, humans often mount an excessive inflammatory response that contributes to chronic lesions, whereas mice, despite sharing several immunological features with humans, develop a more functional response that allows faster lesion resolution.
Our temporal analysis in mice captured, for the first time in this model, the earliest transcriptional responses to infection (2–48 h and Day 14). A peak of immune activation was already evident at 2 h (Module M3), involving antigen presentation, innate immunity, and leukocyte recruitment (Figures 3, 5, Supplementary Figure 2, Appendix). This was followed by the induction of a wound-healing program enriched in keratin and keratin-associated genes (Module M6), which peaked at 6 h and re-emerged at Day 14, preceding ulcer formation (Figure 5, Supplementary Figure 5, Appendix).
By 48 h, however, most transcriptional modules returned close to baseline levels, resulting in a “silent phase” that persisted until Day 14. PCA, DEG profiles, and GSEA analyses all support the presence of this quiescent window. While the functional consequences of this transcriptional state cannot be directly inferred from the present data, it may reflect a period of limited immune activation that temporally coincides with early parasite persistence, as described in other Leishmania infection models (55). At Day 14, the renewed induction of keratin-associated genes (M6) coincided with downregulation of collagen/EMT-related genes (Module M4), potentially priming the tissue for lesion opening (Figure 5).
By Day 35, when ulcerated lesions were established, epidermal and skin barrier genes emerged as a major point of divergence between species. In mice, structural proteins such as FLG, FLG2, and LOR, along with stress-response genes including Crnn and Rptn, which collectively promote keratinocyte differentiation and cornified envelope formation (56–60) were strongly upregulated (Figure 2). Likewise, late cornified envelope (LCE) members (Lce3b, Lce3c, Lce1f) (61, 62) and desmosome components (Dsc1-2), which contribute to repair during inflammation and maintain epidermal cohesion, showed a general pattern of activation in mice but not in human lesions.
Keratinocyte-associated KRT/KRTP genes were broadly and persistently expressed in mice (Supplementary Figure 5, Appendix), consistent with recent single-cell studies indicating that murine migratory keratinocytes retain proliferative potential, whereas in humans migratory and proliferative keratinocytes form distinct, non-overlapping populations (20).
Taken together, these findings suggest that mice engage an early compensatory repair program that preserves barrier integrity and supports rapid epidermal closure (63, 64). In contrast, humans fail to sustain equivalent epidermal and barrier gene expression, resulting in progressive barrier deterioration, chronic inflammation, and impaired healing (Figure 7).
4.3 Comparison with L. major lesions reveals common mechanisms of parasite survival
A recent single-cell study of L. major lesions in C57BL/6 mice revealed broad transcriptional changes and cellular heterogeneity 28 days post-infection (10). Since this model also develops self-healing, Th1-driven lesions (65), we compared it to our L. braziliensis dataset using GSEA (Figure 3). The L. major Day-28 profile most closely resembled an intermediate stage in the L. braziliensis response, lying between Day 14 and Day 35, and all pathways upregulated in L. major at Day 28 were enriched in L. braziliensis at Day 35.
In contrast, downregulation of Ribosome-related pathways emerged earlier in L. braziliensis, during the pre-ulcerative phase at Day 14. Although the functional consequences of this suppression of protein synthesis remain to be experimentally validated, it may reflect a state of reduced host cellular activity that is permissive for parasite persistence during early lesion development, as suggested in other intracellular infection models (66, 67). In support of this interpretation, most pathways at Day 14 were negatively regulated in both KEGG and Reactome databases, with keratinization and cornified envelope formation being the main exceptions (Figure 3).
Reduced protein synthesis at this stage may also impair antigen presentation, further dampening host responses (68), although direct functional validation will be required to substantiate this possibility. In line with this interpretation, alterations in host protein synthesis machinery have also been reported in independent in vivo models of Leishmania infection. Specifically, a study by Venugopal et al. (10) analyzing L. major lesions in the ears of C57BL/6 mice at Day 28 post-infection described a coordinated reduction in the expression of multiple ribosomal subunits, accompanied by decreased signaling through the eukaryotic initiation factor 2 (EIF2) pathway. While that analysis was limited to a single experimental time point during the ulcerated phase, our longitudinal analysis of L. braziliensis infection in BALB/c mice reveals that the most pronounced inhibition of Ribosome- and Translation-related pathways occurs earlier, at Day 14 post-infection, prior to overt lesion ulceration. At this time point, we observed a broad transcriptional suppression characterized by the downregulation of 88 pathways and the upregulation of only four pathways, consistent with a global reduction in host cellular activity. Together, these findings suggest that suppression of host protein synthesis represents a conserved and temporally regulated feature of Leishmania infection, while underscoring the need for additional functional studies to determine whether this transcriptional state reflects active parasite-driven modulation or a host-mediated stress response.
Overall, these cross-species comparisons highlight conserved molecular signatures associated with parasite survival and lesion progression, reinforcing the value of longitudinal and comparative transcriptomic approaches for dissecting host–parasite interactions.
4.4 Skin draining lymph nodes and lesions cross talk
The dLN is a key secondary lymphoid organ that initiates Th1 responses following antigen drainage from the skin (65). To date, no transcriptomic study has evaluated the temporal crosstalk between lesions and corresponding dLNs.
Our data indicate that immune responses, and ultimately disease “resolution,” occur more rapidly in the dLN than in the skin. By Day 77, no DEGs were detected in the dLN, and at Day 35, downregulated genes predominated, consistent with a shift toward immune contraction (Supplementary Figure 9, Appendix).
Communication between the two tissues is reflected in the coordinated regulation of immunoglobulin expression, peaking at Day 35 in the dLN and at Day 77 in the skin (Supplementary Figure 7, Appendix), in line with increasing B-cell frequencies identified by deconvolution analyses (Supplementary Figures 7A, 11B, Appendix). A concurrent reduction in CD4+ T cells in the dLN at Day 35 likely reflects trafficking of effector cells to the lesion site (Supplementary Figure 11, Appendix).
Some inflammatory markers showed divergent regulation between tissues. For instance, S100a8 and S100a9, along with phagocyte and neutrophil markers, were downregulated in the dLN but upregulated in lesions (Supplementary Figure 10D, Appendix). Since these alarmins are important early inflammatory mediators (69) and their absence impairs anti-Leishmania responses (70), their downregulation in dLN beginning at Day 14 may reflect either parasite-driven modulation or changes in immune cell composition.
4.5 Role of immunoglobulins in lesion resolution in the murine model
In human CL caused by L. braziliensis, lesions containing parasite transcripts (PTPos) show elevated B-cell and immunoglobulin expression together with immune inhibitory molecules, suggesting that humoral responses may contribute to parasite persistence (6).
In mice, however, immunoglobulin expression peaks later, at Day 77, when lesions are clinically “resolved” and no active parasite transcription is detected (Supplementary Figure 1C, Appendix). This coincides with strong expression of B cell markers (CD79A, CD19, MS4A1) and activation of a dedicated B-cell module (M5) (Figure 2, Supplementary Figure 7, Appendix). By contrast, during the ulcerated phase (Day 35), lesions show intense inflammation and active parasite transcription but low immunoglobulin transcript levels (Supplementary Figures 1C and 7D, Appendix).
Together, these findings suggest that in humans, immunoglobulin expression overlaps with parasite persistence, likely reflecting less efficient effector mechanisms, whereas in mice, parasites are cleared before the humoral response peaks (Figure 7). Thus, high immunoglobulins levels do not appear to impair lesion resolution in murine L. braziliensis model.
4.6 ncRNAs and epigenetic regulators as host-directed targets in cutaneous leishmaniasis
A notable finding was the high proportion of downregulated ncRNAs at Day 35, accounting for approximately 26% of all DEGs. Robust ncRNA modulation has also been reported in in vitro infection models (71). Several miRNAs detected here, including members of the let-7 family and miR-155, were previously described in vitro, and our data now confirm their regulation in vivo. In addition, we identified a set of miRNAs, such as miR-23a and miR-99a, as novel candidates linked to Leishmania infection (Figure 6A).
Among miRNAs consistently regulated in both human and murine ulcerated lesions, therapeutic prioritization must balance the resolution of immunopathology with the preservation of effective microbicidal responses. Within a host-directed therapy framework, we emphasize modulation of excessive inflammation, particularly as these interventions are envisioned as adjuncts to antiparasitic chemotherapy, thereby reducing the risk of uncontrolled Leishmania replication. Importantly, the therapeutic impact of HDT modulation is likely to be time-sensitive, with efficacy relying on precise intervention windows rather than on constitutive inhibition or activation.
miR-155 and miR-142, both upregulated and associated with pro-inflammatory immune programs, emerged as compelling but complex candidates for antagomir-based inhibition. miR-155, a key regulator of effector CD8+ T cell responses, amplifies cytotoxicity and inflammatory signaling (72) and may contribute to tissue damage in chronic L. braziliensis lesions. Its inhibition could therefore theoretically attenuate CD8+ T cell driven immunopathology. However, in L. major dermal infection model, miR-155 promotes disease susceptibility by favoring Th2 polarization and impairing DC function, while being dispensable for Th1 induction and IFN-γ production by NK and CD8+ T cells (45). Whether miR-155 inhibition in L. braziliensis would primarily reduce pathogenic inflammation or compromise host defense remains unresolved and requires preclinical evaluation. Similarly, miR-142 has broad immunological roles that warrant careful consideration in therapeutic modulation. It regulates T lymphocytes, NK cells, and dendritic cells, and mir142 deficiency results in impaired hematopoiesis. miR-142 also supports neutrophil chemotaxis and effective bacterial clearance during skin infection and wound repair (73). Thus, although its inhibition may limit inflammatory cell recruitment and tissue damage, it may also impair protective immune mechanisms, underscoring the need for careful preclinical evaluation in leishmanial disease.
Interpretation of downregulated miRNAs likewise requires caution, as several are transiently suppressed during normal wound repair. The downregulation of miR-99a, miR-100, and miR-23b represents permissive, pro-healing events during the physiological course of tissue repair (74, 75). Accordingly, further suppression of these miRNAs beyond their physiological downregulation may therefore be counterproductive. In the case of miR-205, early studies reported increased expression during the granulation tissue formation phase of cutaneous wound healing (76). However, subsequent in situ hybridization analyses revealed strong miR-205 expression in the hyperplastic epidermis surrounding the wound, but reduced levels in the leading migrating keratinocytes of the epithelial tongue (77). While it is tempting to speculate that mimic miR-205 restoration could promote wound repair, its spatially restricted expression suggests cell type–specific functions that cannot be inferred from bulk transcriptomic data alone. Approaches such as in situ hybridization or spatial transcriptomics will therefore be essential to accurately define its cellular sources and therapeutic relevance in chronic leishmanial lesions. Together, these miRs candidates highlight mechanistically relevant pathways, but also underscore the challenges of translating them into feasible targets for host-directed therapy in cutaneous leishmaniasis (Figure 7).
Epigenetic regulators also emerged as conserved inflammatory mediators. Sp140 and Kdm6b were upregulated in both murine and human lesions, indicating conserved roles in Leishmania-induced inflammation (Figure 6B). SP140, a bromodomain-containing epigenetic reader, promotes pro-inflammatory transcription and macrophage activation (78, 79). Its inhibitor, GSK761, reduces inflammatory gene expression in macrophages and shows positive effects in Crohn’s disease (80). KDM6B, a lysine demethylase involved in macrophage and DC activation (81, 82), is inhibited by GSK-J4, which reduces cytokine production, including TNF-α, and improves outcomes in models of rheumatoid arthritis, encephalomyelitis and sepsis (83–85). Because epigenetic regulators exert broad transcriptional control, their systemic targeting raises legitimate concerns regarding toxicity, pleiotropy, and cell-type specificity. In the context of cutaneous leishmaniasis, however, these limitations may be mitigated through localized, lesional delivery rather than systemic administration, enabling transient and spatially confined modulation of pathogenic inflammation. Taken together, these results support SP140 and KDM6B as conserved, host-directed targets that could be combined with leishmanicidal drugs to restrain inflammation and improve lesion resolution (Figure 7).
Some limitations should be considered when interpreting our findings. Differences in sampling strategy, with human punch biopsies collected at the lesion edge versus whole ear pinna samples in mice, may influence the relative representation of immune and stromal cell populations. In addition, bulk RNA-seq averages transcriptional signals across heterogeneous cell populations, potentially masking cell type–specific or spatially restricted responses that would be better resolved using single-cell or spatial approaches. Reference-based immune cell deconvolution further relies on single-cell datasets derived from acute cutaneous wounds, which may not fully capture the cellular complexity or activation states present in chronic Leishmania lesions; accordingly, inferred macrophage populations should be interpreted as relative functional programs rather than fixed polarization states. The modest sample size (n = 3 per time point), together with the lack of orthogonal transcript or protein-level validation, may also limit the detection and confirmation of subtle effects, particularly at early stages. Consequentely, key genes and pathways identified here should be regarded as prioritized candidates for future functional validation.
Finally, as RNA extraction was not optimized for small RNA recovery, the detected miRNA profiles likely represent only a subset of the broader miRNA response and should be considered candidate regulators rather than an exhaustive catalog. Despite these constraints, the strong concordance of transcriptional programs observed in ulcerated lesions across species supports the robustness of the core findings and provides a solid framework for future studies. These include analyses of lesions caused by other Leishmania species (e.g., L. amazonensis) and the application of single-cell or spatial transcriptomics to resolve lesion architecture and host–parasite interactions in greater detail. Recent applications of spatial transcriptomics to human cutaneous leishmaniasis biopsies (86, 87) represent an important advance, although high-resolution spatial analyses remain limited. Comparative high-resolution spatial transcriptomics of human and mouse lesions would establish a powerful framework to dissect lesion progression and resolution by integrating parasite distribution, cellular composition and molecular signatures.
In summary, our temporal transcriptomics analysis of L. braziliensis infection in mice reveals that this model faithfully reproduces the core inflammatory and immunopathological features of human cutaneous leishmaniasis, while mounting a more efficient effector response that accelerates lesion resolution. This response is characterized by the concurrent expression of Il12 and Il4, Arg1 and Nos2, and by higher levels of Gzma relative to Gzmb. Beyond inflammation, we show that the differential regulation of epidermal and barrier-related genes underlies the divergent kinetics of wound healing between species. Finally, cross-species comparison of ulcerated lesions identified a set of conserved non-coding RNAs and epigenetic regulators—such as miR-155, miR-142, SP140, and KDM6B—that emerge as promising host-directed targets for therapeutic intervention and warrant further functional validation in pre-clinical models.
Statements
Data availability statement
The datasets generated for this study can be found in the NCBI’s Sequence Read Archive—SRA (https://www.ncbi.nlm.nih.gov/sra/) repository under BioProject ID PRJNA974386. All other relevant data are within the article and its Supporting Information files.
Ethics statement
The animal study was approved by Ethics Committee on the Use of Animals (Comissão de Ética no Uso de Animais – CEUA), FIOCRUZ - Instituto Gonçalo Moniz, Salvador, Bahia-Brazil. The study was conducted in accordance with the local legislation and institutional requirements.
Author contributions
JL-S: Data curation, Methodology, Formal analysis, Writing – original draft, Investigation, Writing – review & editing. CO: Writing – review & editing, Methodology, Investigation. AN: Methodology, Writing – review & editing, Formal analysis, Software. BR: Visualization, Writing – review & editing, Methodology, Data curation, Software. RF: Formal analysis, Writing – original draft, Writing – review & editing, Investigation. JS: Methodology, Software, Writing – review & editing, Data curation, Formal analysis. SS: Writing – review & editing, Methodology, Investigation. SN: Visualization, Writing – review & editing. VM-J: Investigation, Methodology, Writing – review & editing. AB: Writing – review & editing, Funding acquisition, Resources. NT: Writing – review & editing, Data curation. PR: Conceptualization, Writing – review & editing, Investigation, Writing – original draft, Data curation, Funding acquisition, Formal analysis, Visualization. RK: Data curation, Supervision, Formal analysis, Conceptualization, Writing – original draft, Funding acquisition, Writing – review & editing. LF: Funding acquisition, Project administration, Formal analysis, Writing – review & editing, Data curation, Supervision, Writing – original draft, Investigation, Conceptualization, Methodology.
Funding
The author(s) declared that financial support was received for this work and/or its publication. This work was supported by grants from Fundação Oswaldo Cruz (Fiocruz) to LP (PIAP/IGM/FIOCRUZ BA/001/2017 and Inova Fiocruz Terapias Avançadas VPPIS-004-FIO-22) and from Fundação de Amparo a Pesquisa do Estado da Bahia (FAPESB) to AB (PNX0007/2014). This study was also financed in part by the Coordenação de Aperfeiçoamento de Pessoal de Nível Superior (CAPES), Brazil -Finance Code 001. Fellowships were awarded by CAPES to JL-S, by FAPESB to BR, and by Conselho Nacional de Desenvolvimento Científico e Tecnológico (CNPq) to CO, AN, JS, SN, AB, NT, PR and ARK.
Acknowledgments
We are grateful to Dra. Anna Christina de Matos Salim and Dr. Jeronimo (IRR-Fiocruz) for their assistance with library preparation and sequencing of lesion samples, and to Dr. Luiz Lehmann Coutinho (Centro de Genômica Funcional, ESALQ – USP) for his support with library preparation and sequencing of draining lymph node samples. We thank Drs. Artur Trancoso Lopo Queiroz and Theolis Costa Barbosa Bessa for their institutional support and contributions during the initial stages of proposal development and submission. Finally, we would like to acknowledge NEGP-Fiocruz-BA for support in resource management. Figures 1, 4 and 7 were created using BioRender.com.
Conflict of interest
The author(s) declared that this work was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.
Generative AI statement
The author(s) declared that generative AI was used in the creation of this manuscript. Generative AI assistance (ChatGPT, OpenAI) was used solely for grammar and wording revisions.
Any alternative text (alt text) provided alongside figures in this article has been generated by Frontiers with the support of artificial intelligence and reasonable efforts have been made to ensure accuracy, including review by the authors wherever possible. If you identify any issues, please contact us.
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/fimmu.2026.1752217/full#supplementary-material
References
1
KayePScottP. Leishmaniasis: complexity at the host-pathogen interface. Nat Rev Microbiol. (2011) 9:604–15. doi: 10.1038/nrmicro2608. PMID:
2
Cutaneous leishmaniasis | DNDi (2020). Available online at: https://dndi.org/diseases/cutaneous-leishmaniasis/ (Accessed February 3, 2026).
3
AlvarJVélezIDBernCHerreroMDesjeuxPCanoJet al. Leishmaniasis worldwide and global estimates of its incidence. PloS One. (2012) 7:e35671. doi: 10.1371/journal.pone.0035671. PMID:
4
Ministério da Saúde. Manual de vigilância da Leishmaniose tegumentar. Brasília: Ministério da Saúde (2017). 189. Available online at: https://bvsms.saude.gov.br/bvs/publicacoes/manual_vigilancia_leishmaniose_tegumentar.pdf.
5
NovaisFOCarvalhoLPPassosSRoosDSCarvalhoEMScottPet al. Genomic profiling of human Leishmania Braziliensis lesions identifies transcriptional modules associated with cutaneous immunopathology. J Invest Dermatol. (2015) 135:94–101. doi: 10.1038/jid.2014.305. PMID:
6
ChristensenSMDillonLALCarvalhoLPPassosSNovaisFOHughittVKet al. Meta-transcriptome profiling of the human-Leishmania Braziliensis cutaneous lesion. PloS NeglTrop Dis. (2016) 10:e0004992. doi: 10.1371/journal.pntd.0004992. PMID:
7
AmorimCFNovaisFONguyenBTMisicAMCarvalhoLPCarvalhoEMet al. Variable gene expression and parasite load predict treatment outcome in cutaneous leishmaniasis. Sci Transl Med. (2019) 11:eaax4204. doi: 10.1126/scitranslmed.aax4204. PMID:
8
Farias AmorimCO. NovaisFNguyenBTNascimentoMTLagoJLagoASet al. Localized skin inflammation during cutaneous leishmaniasis drives a chronic, systemic IFN-γ signature. PloS NeglTrop Dis. (2021) 15:e0009321. doi: 10.1371/journal.pntd.0009321. PMID:
9
Farias AmorimCLovinsVMSinghTPNovaisFOHarrisJCLagoASet al. Multiomic profiling of cutaneous leishmaniasis infections reveals microbiota-driven mechanisms underlying disease severity. Sci Transl Med. (2023) 15. doi: 10.1126/scitranslmed.adh1469. PMID:
10
VenugopalGBirdJTWashamCLRoysHBowlinAByrumSDet al. In vivo transcriptional analysis of mice infected with Leishmania major unveils cellular heterogeneity and altered transcriptomic profiling at single-cell resolution. PloS NeglTrop Dis. (2022) 16:e0010518. doi: 10.1371/journal.pntd.0010518. PMID:
11
VenugopalGBirdJTRoysHBowlinAFryLByrumSDet al. Both the infection status and inflammatory microenvironment induce transcriptional remodeling in macrophages in murine leishmanial lesions. J Parasitol. (2023) 109:200–10. doi: 10.1645/22-94. PMID:
12
LeeSHKangBKamenyevaOFerreiraTRChoKKhillanJSet al. Dermis resident macrophages orchestrate localized ILC2 eosinophil circuitries to promote non-healing cutaneous leishmaniasis. Nat Commun. (2023) 14:7852. doi: 10.1038/s41467-023-43588-2. PMID:
13
De MouraTRNovaisFOOliveiraFClarêncioJNoronhaABarralAet al. Toward a novel experimental model of infection to study American cutaneous leishmaniasis caused by Leishmania Braziliensis. Infect Immun. (2005) 73:5827–34. doi: 10.1128/IAI.73.9.5827-5834.2005. PMID:
14
NovaisFOAmorimCFScottP. Host-directed therapies for cutaneous leishmaniasis. Front Immunol. (2021) 12:660183. doi: 10.3389/fimmu.2021.660183. PMID:
15
LecoeurHPrinaEGutiérrez-SanchezMSpäthGF. Going ballistic: Leishmania nuclear subversion of host cell plasticity. Trends Parasitol. (2022) 38:205–16. doi: 10.1016/j.pt.2021.09.009. PMID:
16
SanzCRMiróGSevaneNReyes-PalomaresADunnerS. Modulation of host immune response during Leishmania infantum natural infection: a whole-transcriptome analysis of the popliteal lymph nodes in dogs. Front Immunol. (2022) 12:794627. doi: 10.3389/fimmu.2021.794627. PMID:
17
NgH-ITuongZKFernandoGJPDepelsenaireACIMeligaSCFrazerIHet al. Microprojection arrays applied to skin generate mechanical stress, induce an inflammatory transcriptome and cell death, and improve vaccine-induced immune responses. NPJ Vaccines. (2019) 4:1–12. doi: 10.1038/s41541-019-0134-4. PMID:
18
ReynoldsGVeghPFletcherJPoynerEFMStephensonEGohIet al. Developmental cell programs are co-opted in inflammatory skin disease. Science. (2021) 371. doi: 10.1126/science.aba6500. PMID:
19
CaiYXiongMXinZLiuCRenJYangXet al. Decoding aging-dependent regenerative decline across tissues at single-cell resolution. Cell Stem Cell. (2023) 30:1674–1691.e8. doi: 10.1016/j.stem.2023.09.014. PMID:
20
LiuZBianXLuoLBjörklundÅKLiLZhangLet al. Spatiotemporal single-cell roadmap of human skin wound healing. Cell Stem Cell. (2025) 32:479–498.e8. doi: 10.1016/j.stem.2024.11.013. PMID:
21
BolgerAMLohseMUsadelB. Trimmomatic: a flexible trimmer for Illumina sequence data. Bioinforma Oxf Engl. (2014) 30:2114–20. doi: 10.1093/bioinformatics/btu170. PMID:
22
FrankishACarbonell-SalaSDiekhansMJungreisILovelandJEMudgeJMet al. GENCODE: reference annotation for the human and mouse genomes in 2023. Nucleic Acids Res. (2023) 51:D942–9. doi: 10.1093/nar/gkac1071. PMID:
23
DobinADavisCASchlesingerFDrenkowJZaleskiCJhaSet al. STAR: ultrafast universal RNA-seq aligner. Bioinforma Oxf Engl. (2013) 29:15–21. doi: 10.1093/bioinformatics/bts635. PMID:
24
la FuenteSGCamachoEPeiró-PastorRRastrojoACarrasco-RamiroFAguadoBet al. Complete and de novo assembly of the Leishmania Braziliensis (M2904) genome. Mem Inst Oswaldo Cruz. (2018) 114:e180438. doi: 10.1590/0074-02760180438. PMID:
25
RobinsonMDOshlackA. A scaling normalization method for differential expression analysis of RNA-seq data. Genome Biol. (2010) 11:R25. doi: 10.1186/gb-2010-11-3-r25. PMID:
26
RobinsonMDMcCarthyDJSmythGK. edgeR: a Bioconductor package for differential expression analysis of digital gene expression data. Bioinforma Oxf Engl. (2010) 26:139–40. doi: 10.1093/bioinformatics/btp616. PMID:
27
SubramanianATamayoPMoothaVKMukherjeeSEbertBLGilletteMAet al. Gene set enrichment analysis: a knowledge-based approach for interpreting genome-wide expression profiles. Proc Natl Acad Sci. (2005) 102:15545–50. doi: 10.1073/pnas.0506580102. PMID:
28
RamosPIPCristalJRKhouriRBoaventuraVAzevedoLGCorreiaTCet al. Selective suppression of cellular immunity and increased cytotoxicity in skin lesions of disseminated leishmaniasis uncovered by transcriptome-wide analysis. J Invest Dermatol. (2021) 141:2542–2546.e5. doi: 10.1016/j.jid.2021.03.017. PMID:
29
KanehisaM. KEGG: kyoto encyclopedia of genes and genomes. Nucleic Acids Res. (2000) 28:27–30. doi: 10.1093/nar/28.1.27. PMID:
30
MilacicMBeaversDConleyPGongCGillespieMGrissJet al. The reactome pathway knowledgebase 2024. Nucleic Acids Res. (2024) 52:D672–8. doi: 10.1093/nar/gkad1025. PMID:
31
ChenZQuanLHuangAZhaoQYuanYYuanXet al. seq-ImmuCC: cell-centric view of tissue transcriptome measuring cellular compositions of immune microenvironment from mouse RNA-seq data. Front Immunol. (2018) 9:1286. doi: 10.3389/fimmu.2018.01286. PMID:
32
WickhamH. ggplot2: elegant graphics for data analysis. Cham: Springer international publishing (2016). p. 1.
33
DaiYHuRLiuAChoKSManuelAMLiXet al. WebCSEA: web-based cell-type-specific enrichment analysis of genes. Nucleic Acids Res. (2022) 50:W782–90. doi: 10.1093/nar/gkac392. PMID:
34
RussoPSTFerreiraGRCardozoLEBürgerMCArias-CarrascoRMaruyamaSRet al. CEMiTool: a Bioconductor package for performing comprehensive modular co-expression analyses. BMC Bioinf. (2018) 19:56. doi: 10.1186/s12859-018-2053-1. PMID:
35
BreuerKForoushaniAKLairdMRChenCSribnaiaALoRet al. InnateDB: systems biology of innate immunity and beyond--recent updates and continuing curation. Nucleic Acids Res. (2013) 41:D1228–33. doi: 10.1093/nar/gks1147. PMID:
36
DasSMcClainCJRaiSN. Fifteen years of gene set analysis for high-throughput genomic data: a review of statistical approaches and future challenges. Entropy. (2020) 22:427. doi: 10.3390/e22040427. PMID:
37
Van DamSVõsaUVan Der GraafAFrankeLDe MagalhãesJP. Gene co-expression analysis for functional classification and gene–disease predictions. Brief Bioinform. (2017) 19(4):575–92. doi: 10.1093/bib/bbw139. PMID:
38
FariasLP. Interactive HTML figure for Temporal transcriptional dynamics in cutaneous leishmaniasis reveal novel targets for therapeutic interventions in a dermal mouse model (2025). Available online at: https://lpfarias.github.io/cemitools-Leish-lesion/ (Accessed February 3, 2026).
39
FariasLP. Interactive HTML figure for Temporal transcriptional dynamics in cutaneous leishmaniasis reveal novel targets for therapeutic interventions in a dermal mouse model (2025). Available online at: https://lpfarias.github.io/cemitools-Leish-sdLN/ (Accessed February 3, 2026).
40
KumarAVijaykumarSDikhitMRAbhishekKMukherjeeRSenAet al. Differential regulation of miRNA profiles of human cells experimentally infected by Leishmania donovani isolated from Indian visceral leishmaniasis and post-kala-azar dermal leishmaniasis. Front Microbiol. (2020) 11:1716. doi: 10.3389/fmicb.2020.01716. PMID:
41
Linhares-LacerdaLTemerozoJRRibeiro-AlvesMAzevedoEPMojoliANascimentoMTCet al. Neutrophil extracellular trap-enriched supernatants carry microRNAs able to modulate TNF-α production by macrophages. Sci Rep. (2020) 10:2715. doi: 10.1038/s41598-020-59486-2. PMID:
42
MendonçaLSOSantosJMKanetoCMde CarvalhoLDLima-SantosJAugustoDGet al. Characterization of serum cytokines and circulating microRNAs that are predicted to regulate inflammasome genes in cutaneous leishmaniasis patients. Exp Parasitol. (2020) 210:107846. doi: 10.1016/j.exppara.2020.107846. PMID:
43
MuxelSMAcuñaSMAokiJIZampieriRAFloeter-WinterLM. Toll-like receptor and miRNA-let-7e expression alter the inflammatory response in Leishmania amazonensis-infected macrophages. Front Immunol. (2018) 9:2792. doi: 10.3389/fimmu.2018.02792. PMID:
44
SilvaSCSilvaDFAlmeidaTCPerasoliFBDa SilvaATPDa SilvaGNet al. Behavior of two Leishmania infantum strains—evaluation of susceptibility to antimonials and expression of microRNAs in experimentally infected J774 macrophages and in BALB/c mice. Parasitol Res. (2018) 117:2881–93. doi: 10.1007/s00436-018-5979-3. PMID:
45
VarikutiSVermaCNatarajanGOghumuSSatoskarAR. MicroRNA155 plays a critical role in the pathogenesis of cutaneous Leishmania major infection by promoting a Th2 response and attenuating dendritic cell activity. Am J Pathol. (2021) 191:809–16. doi: 10.1016/j.ajpath.2021.01.012. PMID:
46
BiswasSRaoCM. Epigenetic tools (the writers, the readers and the erasers) and their implications in cancer therapy. Eur J Pharmacol. (2018) 837:8–24. doi: 10.1016/j.ejphar.2018.08.021. PMID:
47
SacksDNoben-TrauthN. The immunology of susceptibility and resistance to Leishmania major in mice. Nat Rev Immunol. (2002) 2:845–58. doi: 10.1038/nri933. PMID:
48
BrodskynCde OliveiraCIBarralABarral-NettoM. Vaccines in leishmaniasis: advances in the last five years. Expert Rev Vaccines. (2003) 2:705–17. doi: 10.1586/14760584.2.5.705. PMID:
49
ScorzaBMCarvalhoEMWilsonME. Cutaneous manifestations of human and murine leishmaniasis. Int J Mol Sci. (2017) 18:1296. doi: 10.3390/ijms18061296. PMID:
50
CardosoTMMaChado>ÁCostaDLCarvalhoLPQueirozAMaChadoPet al. Protective and pathological functions of CD8+ T cells in Leishmania Braziliensis infection. Infect Immun. (2015) 83:898–906. doi: 10.1128/IAI.02404-14. PMID:
51
SantosCSBoaventuraVRibeiro CardosoCTavaresNLordeloMJNoronhaAet al. CD8(+) granzyme B(+)-mediated tissue injury vs. CD4(+)IFNγ(+)-mediated parasite killing in human cutaneous leishmaniasis. J Invest Dermatol. (2013) 133:1533–40. doi: 10.1038/jid.2013.4. PMID:
52
NovaisFOCarvalhoAMClarkMLCarvalhoLPBeitingDPBrodskyIEet al. CD8+ T cell cytotoxicity mediates pathology in the skin by inflammasome activation and IL-1β production. PloS Pathog. (2017) 13:e1006196. doi: 10.1371/journal.ppat.1006196. PMID:
53
FowlerEAFarias AmorimCMostacadaKYanAAmorim SacramentoLStancoRAet al. Neutrophil-mediated hypoxia drives pathogenic CD8+ T cell responses in cutaneous leishmaniasis. J Clin Invest. (2024) 134:e177992. doi: 10.1172/JCI177992. PMID:
54
PeñaOAMartinP. Cellular and molecular mechanisms of skin wound healing. Nat Rev Mol Cell Biol. (2024) 25:599–616. doi: 10.1038/s41580-024-00715-1. PMID:
55
BelkaidYMendezSLiraRKadambiNMilonGSacksD. A natural model of Leishmania major infection reveals a prolonged “silent” phase of parasite amplification in the skin before the onset of lesion formation and immunity. J Immunol. (2000) 165:969–77. doi: 10.4049/jimmunol.165.2.969. PMID:
56
SandilandsASutherlandCIrvineADMcLeanWHI. Filaggrin in the frontline: role in skin barrier function and disease. J Cell Sci. (2009) 122:1285–94. doi: 10.1242/jcs.033969. PMID:
57
CandiESchmidtRMelinoG. The cornified envelope: a model of cell death in the skin. Nat Rev Mol Cell Biol. (2005) 6:328–40. doi: 10.1038/nrm1619. PMID:
58
MehrelTHohlDRothnagelJALongleyMABundmanDChengCet al. Identification of a major keratinocyte cell envelope protein, loricrin. Cell. (1990) 61:1103–12. doi: 10.1016/0092-8674(90)90073-N. PMID:
59
ContzlerRFavreBHuberMHohlD. Cornulin, a new member of the “fused gene” family, is expressed during epidermal differentiation. J Invest Dermatol. (2005) 124:990–7. doi: 10.1111/j.0022-202X.2005.23694.x. PMID:
60
KriegPSchupplerMKoestersRMinchevaALichterPMarksF. Repetin (Rptn), a new member of the “fused gene” subgroup within the S100 gene family encoding a murine epidermal differentiation protein. Genomics. (1997) 43:339–48. doi: 10.1006/geno.1997.4818. PMID:
61
JacksonBTilliCMLJHardmanMJAvilionAAMacLeodMCAshcroftGSet al. Late cornified envelope family in differentiating epithelia—response to calcium and ultraviolet irradiation. J Invest Dermatol. (2005) 124:1062–70. doi: 10.1111/j.0022-202X.2005.23699.x. PMID:
62
de CidRRiveira-MunozEZeeuwenPLJMRobargeJLiaoWDannhauserENet al. Deletion of the late cornified envelope LCE3B and LCE3C genes as a susceptibility factor for psoriasis. Nat Genet. (2009) 41:211–5. doi: 10.1038/ng.313. PMID:
63
GreenKJSimpsonCL. Desmosomes: New perspectives on a classic. J Invest Dermatol. (2007) 127:2499–515. doi: 10.1038/sj.jid.5701015. PMID:
64
MüllerLHatzfeldMKeilR. Desmosomes as signaling hubs in the regulation of cell behavior. Front Cell Dev Biol. (2021) 9:745670. doi: 10.3389/fcell.2021.745670. PMID:
65
ScottPNovaisFO. Cutaneous leishmaniasis: immune responses in protection and pathogenesis. Nat Rev Immunol. (2016) 16:581–92. doi: 10.1038/nri.2016.72. PMID:
66
ToribioRVentosoI. Inhibition of host translation by virus infection in vivo. Proc Natl Acad Sci. (2010) 107:9837–42. doi: 10.1073/pnas.1004110107. PMID:
67
ShresthaNBahnanWWileyDJBarberGFieldsKASchesserK. Eukaryotic initiation factor 2 (eIF2) signaling regulates proinflammatory cytokine expression and bacterial invasion. J Biol Chem. (2012) 287:28738–44. doi: 10.1074/jbc.M112.375915. PMID:
68
ArgüelloRJReverendoMGattiEPierreP. Regulation of protein synthesis and autophagy in activated dendritic cells: implications for antigen processing and presentation. Immunol Rev. (2016) 272:28–38. doi: 10.1111/imr.12427. PMID:
69
VoglTEisenblätterMVöllerTZenkerSHermannSVan LentPet al. Alarmin S100A8/S100A9 as a biomarker for molecular imaging of local inflammatory activity. Nat Commun. (2014) 5:4593. doi: 10.1038/ncomms5593. PMID:
70
ContrerasIShioMTCesaroATessierPAOlivierM. Impact of neutrophil-secreted myeloid related proteins 8 and 14 (MRP 8/14) on leishmaniasis progression. PloS NeglTrop Dis. (2013) 7:e2461. doi: 10.1371/journal.pntd.0002461. PMID:
71
FernandesJCRGonçalvesANAFloeter-WinterLMNakayaHIMuxelSM. Comparative transcriptomic analysis of long noncoding RNAs in Leishmania-infected human macrophages. Front Genet. (2023) 13:1051568. doi: 10.3389/fgene.2022.1051568. PMID:
72
DuddaJCSalaunBJiYPalmerDCMonnotGCMerckEet al. MicroRNA-155 is required for effector CD8+ T cell responses to virus infection and cancer. Immunity. (2013) 38:742–53. doi: 10.1016/j.immuni.2012.12.006. PMID:
73
OlsonWJDerudderE. The miR-142 miRNAs: Shaping the naïve immune system. Immunol Lett. (2023) 261:37–46. doi: 10.1016/j.imlet.2023.07.005. PMID:
74
JinYTymenSDChenDFangZJZhaoYDragasDet al. MicroRNA-99 family targets AKT/mTOR signaling pathway in dermal wound healing. PloS One. (2013) 8:e64434. doi: 10.1371/journal.pone.0064434. PMID:
75
ZhangXYangJZhaoJZhangPHuangX. MicroRNA-23b inhibits the proliferation and migration of heat-denatured fibroblasts by targeting Smad3. PloS One. (2015) 10:e0131867. doi: 10.1371/journal.pone.0131867. PMID:
76
WangTFengYSunHZhangLHaoLShiCet al. miR-21 regulates skin wound healing by targeting multiple aspects of the healing process. Am J Pathol. (2012) 181:1911–20. doi: 10.1016/j.ajpath.2012.08.022. PMID:
77
WangTZhaoNLongSGeLWangASunHet al. Downregulation of miR-205 in migrating epithelial tongue facilitates skin wound re-epithelialization by derepressing ITGA5. Biochim Biophys Acta BBA - Mol Basis Dis. (2016) 1862:1443–52. doi: 10.1016/j.bbadis.2016.05.004. PMID:
78
MehtaSCronkiteDABasavappaMSaundersTLAdiliaghdamFAmatullahHet al. Maintenance of macrophage transcriptional programs and intestinal homeostasis by epigenetic reader SP140. Sci Immunol. (2017) 2:eaag3160. doi: 10.1126/sciimmunol.aag3160. PMID:
79
FraschillaIJeffreyKL. The speckled protein (SP) family: Immunity’s chromatin readers. Trends Immunol. (2020) 41:572–85. doi: 10.1016/j.it.2020.04.007. PMID:
80
GhiboubMKosterJCraggsPDLi YimAYFShillingsAHutchinsonSet al. Modulation of macrophage inflammatory function through selective inhibition of the epigenetic reader protein SP140. BMC Biol. (2022) 20:182. doi: 10.1186/s12915-022-01380-6. PMID:
81
De SantaFTotaroMGProsperiniENotarbartoloSTestaGNatoliG. The histone H3 lysine-27 demethylase Jmjd3 links inflammation to inhibition of polycomb-mediated gene silencing. Cell. (2007) 130:1083–94. doi: 10.1016/j.cell.2007.08.019. PMID:
82
DasNDJungKHChoiMRYoonHSKimSHChaiYG. Gene networking and inflammatory pathway analysis in a JMJD3 knockdown human monocytic cell line. Cell Biochem Funct. (2012) 30:224–32. doi: 10.1002/cbf.1839. PMID:
83
KruidenierLChungCChengZLiddleJCheKJobertyGet al. A selective jumonji H3K27 demethylase inhibitor modulates the proinflammatory macrophage response. Nature. (2012) 488:404–12. doi: 10.1038/nature11262. PMID:
84
DoñasCCarrascoMFritzMPradoCTejónGOsorio-BarriosFet al. The histone demethylase inhibitor GSK-J4 limits inflammation through the induction of a tolerogenic phenotype on DCs. J Autoimmun. (2016) 75:105–17. doi: 10.1016/j.jaut.2016.07.011. PMID:
85
PanYWangJXueYZhaoJLiDZhangSet al. GSKJ4 protects mice against early sepsis via reducing proinflammatory factors and up-regulating MiR-146a. Front Immunol. (2018) 9:2272. doi: 10.3389/fimmu.2018.02272. PMID:
86
ParkashVAshwinHDeySSadlovaJVojtkovaBVan BocxlaerKet al. Safety and reactogenicity of a controlled human infection model of sand fly-transmitted cutaneous leishmaniasis. Nat Med. (2024) 30:3150–62. doi: 10.1038/s41591-024-03146-9. PMID:
87
DeyNSDeySBrownNSenarathneSCampos ReisLSenguptaRet al. IL-32–producing CD8+ memory T cells define immunoregulatory niches in human cutaneous leishmaniasis. J Clin Invest. (2025) 135. doi: 10.1172/jci182040. PMID:
Summary
Keywords
BALB/c mouse model, cutaneous leishmaniasis, immune response, Leishmania braziliensis, microRNAs (miRNAs), wound healing, transcriptomics
Citation
Lobo-Silva J, Orge C, Neto APS, Ribeiro BV, Fávaro RD, Silva JK, Santos SPO, Nunes S, Moitinho-Junior VS, Barral A, Machado NT, Ramos PIP, Khouri R and Farias LP (2026) Temporal transcriptional dynamics in cutaneous leishmaniasis reveal novel targets for therapeutic interventions in a dermal mouse model. Front. Immunol. 17:1752217. doi: 10.3389/fimmu.2026.1752217
Received
22 November 2025
Revised
08 February 2026
Accepted
09 March 2026
Published
12 May 2026
Volume
17 - 2026
Edited by
Debanjan Mukhopadhyay, Presidency University, India
Reviewed by
Shima Hadifar, Pasteur Institute of Iran (PII), Iran
Sajad Rashidi, Khomein University of Medical Sciences, Iran
Updates
Copyright
© 2026 Lobo-Silva, Orge, Neto, Ribeiro, Fávaro, Silva, Santos, Nunes, Moitinho-Junior, Barral, Machado, Ramos, Khouri and Farias.
This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.
*Correspondence: Leonardo Paiva Farias, leonardo.farias@fiocruz.br
Disclaimer
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.