In vivo transcriptional analysis of mice infected with Leishmania major unveils cellular heterogeneity and altered transcriptomic profiling at single-cell resolution

Leishmania parasites cause cutaneous leishmaniasis (CL), a disease characterized by disfiguring, ulcerative skin lesions. Both parasite and host gene expression following infection with various Leishmania species has been investigated in vitro, but global transcriptional analysis following L. major infection in vivo is lacking. Thus, we conducted a comprehensive transcriptomic profiling study combining bulk RNA sequencing (RNA-Seq) and single-cell RNA sequencing (scRNA-Seq) to identify global changes in gene expression in vivo following L. major infection. Bulk RNA-Seq analysis revealed that host immune response pathways like the antigen processing and presentation pathway were significantly enriched amongst differentially expressed genes (DEGs) upon infection, while ribosomal pathways were significantly downregulated in infected mice compared to naive controls. scRNA-Seq analyses revealed cellular heterogeneity including distinct resident and recruited cell types in the skin following murine L. major infection. Within the individual immune cell types, several DEGs indicative of many interferon induced GTPases and antigen presentation molecules were significantly enhanced in the infected ears including macrophages, resident macrophages, and inflammatory monocytes. Ingenuity Pathway Analysis of scRNA-Seq data indicated the antigen presentation pathway was increased with infection, while EIF2 signaling is the top downregulated pathway followed by eIF4/p70S6k and mTOR signaling in multiple cell types including macrophages, blood and lymphatic endothelial cells. Altogether, this transcriptomic profile highlights known recruitment of myeloid cells to lesions and recognizes a potential role for EIF2 signaling in murine L. major infection in vivo.


Introduction
Leishmaniasis is a multifaceted disease caused by different species of obligate intracellular protozoan parasites of the genus Leishmania, belonging to the Trypanosomatid family. Depending on the complex interaction between the species and the host immune system, the disease can vary in severity resulting in a wide spectrum of clinical outcomes that have been classified into the following categories: cutaneous leishmaniasis (CL), mucocutaneous leishmaniasis (MCL), disseminated leishmaniasis or diffuse CL (DCL) where lesions remain localized to skin or mucosal surfaces, and a life-threatening condition called visceral leishmaniasis (VL) where parasites migrate to the internal organs like the liver, spleen and bone marrow [1]. As per the World Health Organization (WHO), leishmaniasis is a major public health problem due to its annual incidence of 1.7 million new cases worldwide [1,2].
Of the species belonging to subgenus Leishmania, L. major is an important etiological agent of CL and possesses clinical and epidemiological importance, especially in parts of Asia, the Middle East, Northern Africa, and Southern Europe [3]. Although CL is not fatal and considered to be a self-healing disease, the development of nodules or papules followed by ulcerations at the site of infection is the hallmark of the disease; importantly, both parasite replication and the host immune response can contribute to the disease pathology [3]. CL produces permanently disfiguring skin lesions that are associated with massive immune cell recruitment, across the blood vascular endothelium, and into the skin where the parasite resides [4][5][6][7][8]. In addition to parasite control, expansion of the lymphatics is required for lesion resolution [9,10]. While CD4 + Th1 cells producing IFNγ are required to activate macrophages to kill parasites via NO, the exacerbated activation and sustained recruitment of immune cells including neutrophils, NK cells, Ly6C + inflammatory monocytes, and CD4 + and CD8 + T lymphocytes induces a chronic inflammatory response; this chronic inflammation leads to tissue necrosis and skin damage, a feature of non-healing lesions [11][12][13].
Myeloid cells such as monocytes and macrophages are the primary permissive host cells for Leishmania parasites, and elimination of parasites by these cells is critical for host resistance [11,[14][15][16]. While the importance of macrophages as host cells for parasites and effector cells involved in parasite killing has been well-characterized, emerging evidence suggests Ly6C + CCR2 + inflammatory monocytes are a preferential target for parasites and serve as the initial myeloid host cell for early parasite replication; however, monocytes also perform effector functions by producing inducible nitric oxide (iNOS) in an IFNγ-dependent manner to kill parasites [14,[16][17][18][19]. In addition to the parasite taking advantage of the Th1-mediated recruitment of inflammatory monocytes for infection, recruited neutrophils and monocyte-derived dendritic cells at the site of infection can also harbor parasites, suggesting these immune cells serve an important function in host-parasite interplay [16,[20][21][22][23][24][25][26].
Leishmania parasites possess a variety of virulence mechanisms [27,28] and use several strategies to evade the host immune response for intracellular survival including modulating the host immune response by altering T cell responses, impeding antigen display by MHCII, and hindering nitric oxide (NO) production [29,30]. NO not only directly kills parasites, but NO also regulates the breadth of the inflammatory response through quorum sensing and restricts the supply of proliferation-permissive host cells including monocyte-derived phagocytes to the site of infection [26,31,32]. Importantly, Leishmania parasites can escape from oxidative burst and they fail to activate optimal macrophage innate immune responses [33,34]. Despite the robust Th1 response and NO production, parasite growth persists in CD206 + tissue-resident dermal macrophages in non-healing L. major infection suggesting cells in inflamed skin respond to soluble tissue mediators in the microenvironment differently. As a result, we set out to identify changes at the transcriptional level following Leishmania infection within different cell types and especially within the hostile tissue microenvironment.
Until recently, many studies examining the host response relied mostly on microarraybased or serial analysis of gene expression tag approaches. Using these approaches, previous studies compared the host gene expression profile between infection with promastigotes and amastigotes [35][36][37][38][39] or different species of Leishmania [40][41][42]. However, microarray-based approaches have technical limitations such as hybridization and cross hybridization artefacts, dye-based detection issues, certain probes cannot be included on the microarrays, and the inability to detect 5' and 3' UTRs boundaries [43,44]. In recent years, RNA-Sequencing (RNA-Seq) has emerged as a powerful tool to study transcriptional changes and also significantly enhanced our understanding of disease pathogenesis due to its high sensitivity. Moreover, transcriptomic profiling using RNA-Seq following infection with various species of Leishmania has been mostly applied to in vitro experiments and some studies have investigated transcriptional changes in the human or murine host as well as the parasite simultaneously [45][46][47][48][49][50][51][52]. A recent study comparing the gene expression profile of cutaneous lesions from L. braziliensis-infected patients with or without treatment revealed most of the differentially expressed transcripts were correlated with cytotoxicity-related pathways and parasite loads [49]. The same group also revealed a consistent myeloid interferon stimulated gene (ISG) signature in skin lesions from L. braziliensis-infected patients using RNA-Seq [53]. Altogether, these studies revealed that the host immune response upregulates transcripts related to both pro-and anti-inflammatory responses during leishmaniasis [54,55]. Although the key biological mechanisms involved in the pathogenesis of CL have been well characterized using both in vitro and in vivo models, many of the mechanistic studies in vivo were investigated in the context of a single cytokine or only examining one to few cell types at a time, which does not reflect the complex inflammatory environment of leishmanial lesions. Additionally, many studies predominantly focused on the immune cells and neglected stromal cells such as the endothelial cells, so we aimed to perform a study that encompasses all the cell types located in lesions that may participate in pathogen control and/or the pathogenesis of disease. Furthermore, none of the above RNA-Seq studies to date have characterized the global transcriptional reprogramming following L. major infection in vivo, which is the most widely used murine model of CL to study disease pathogenesis and parasite-host interactions. The Leishmania field also lacks a comprehensive murine transcriptomic profile that could identify novel transcriptomic changes and new pathways within individual cell types in vivo that could lead to exciting new avenues for therapeutic interventions. Therefore, our study combines both bulk RNA-Seq and scRNA-Seq to assemble a comprehensive dataset that defines how the host response reprograms individual cell types following L. major infection to determine the molecular mediators contributing to the pathogenesis of disease.

Enriched pathways for differentially expressed genes during L. major infection revealed by bulk RNA-Seq
To characterize the transcriptomic landscape during L. major infection in vivo, we applied RNA-Seq analysis on ear samples obtained from mice that were infected with L. major promastigotes and uninfected naive control ears. In this experimental murine model of L. major infection, lesion volume peaks between 3-5 weeks post-infection (p.i.) (S1 Fig) and then the lesion starts to resolve spontaneously about 8 weeks p.i. after parasites have been controlled by a Th1 immune response; full lesion resolution typically occurs about 12 weeks p.i. For parasite burdens, at 4 weeks p.i, an average of 2-2.5x10 9 parasites are present in lesions and parasite numbers decline over time as the lesion resolves. As a result, the host transcriptome was investigated from ear samples collected from experimentally-infected mice at 4 weeks p.i. and compared to naive controls [56].
To compare the gene expression profiles of infected and naive mice, bulk RNA-Seq was performed on L. major-infected ears and naive control ears. Transcriptional analysis revealed that ears from infected mice and naive controls were distinct from one another, as determined by multidimensional scaling (MDS) plot and DEG analysis. MDS plot shows the positions of each sample, with samples from different experimental groups being well separated, and samples from the same experimental group clustering together ( Fig 1A). Therefore, the distance between samples reveals the distinct pattern of gene expression between the infected and naive animals ( Fig 1A). To investigate transcriptomic signatures associated with infection, we carried out DEG analysis between infected mice and naive controls by comparing the RNA-Seq read counts of the various genes and subsequently applying the cut-off criteria. Transcripts that were either increased or decreased in RNA abundance with infection (logCPM><1) were included in the volcano plot showing transcriptional differences observed between infected and naive ears ( Fig 1B). The gene expression profiles derived from the RNA-Seq data were calculated using the RPKM method [57] and a fold change >2 and p<0.05 were considered statistically significant. Of more than 10,800 genes that were detectable in infected ears, we observed that 211 transcripts were increased in RNA abundance and 34 transcripts were decreased in RNA abundance with infection, while 10,014 genes did not show any significant differences between naive and L. major-infected ears ( Fig 1B). Here, it should be noted that Gene Set Enrichment Analysis (GSEA) was carried out against the KEGG pathways and involved a ranked order list which contains both DEG and fold changes for more than 10800 genes. The complete list of genes detected in our RNA-Seq analysis is listed in S1 Table. GSEA using KEGG pathways revealed a total of 276 enriched pathways which includes pathways involved in both disease conditions and molecular signaling networks. Specifically, the antigen processing and presentation pathway was found to be significantly enriched amongst DEGs, while the ribosomal pathway was significantly downregulated during L. major infection (Fig 1C and 1D). In addition to antigen processing and presentation, we observed many other host immune response pathways upregulated with infection including: cytokinecytokine receptor interaction, phagosome, chemokine signaling, cell-adhesion molecules pathway, NK cell mediated cytotoxicity, leukocyte trans-endothelial migration and Fcγ receptormediated phagocytosis (Fig 1C). Conversely, top 20 downregulated pathways enriched for DEGs in L. major infection were related to ribosomal translation, mineral absorption, ECMreceptor interaction, biosynthesis of unsaturated fatty acids, and steroid biosynthesis (Fig 1D).
The top 10 KEGG pathways for both upregulated and downregulated pathways with the fold change and adjusted 'p' value are listed in Table 1. A number of disease-specific KEGG pathways appeared prominent in the enrichment analysis including Staphylococcus aureus infection, autoimmune thyroid disease, and graft vs. host disease (S2 Table). Importantly,  Table).

scRNA-Seq reveals the cellular heterogeneity and altered transcriptomic profile of individual cell types during murine L. major infection in vivo
The bulk RNA-Seq analysis revealed global changes in the transcriptional profile between infected mice and naive controls following L. major inoculation. To further investigate transcriptomic changes within individual cell types present in leishmanial lesions, scRNA-Seq was performed to provide a deeper understanding of how individual cells function in the tissue microenvironment. Single cells from the ears of infected and naive mice were bar-coded and sequenced using the droplet-based 10X Genomics Chromium platform ( Fig 3A). After quality control assessment and filtering, the datasets were processed using Cell Ranger software. Unbiased hierarchical clustering using the Seurat R package provides single-cell transcriptional profiling with 26,558 cells and displayed the cellular heterogeneity which includes both resident and recruited cell types. Cell populations from 35 distinct cell types were defined using canonical markers from published literature and online databases ( Fig 3B) [58]. The dot plot representing the cell type-specific canonical markers for each cell lineage used to distinguish the 35 distinct clusters is provided ( Fig 4A). Here, it is to be noted that we couldn't separate out the free parasites that are possibly being attached to the various cell types. Amongst the 35 cell types, we identified 16 cell types containing immune cells. Feature plots show the expression of cell type-specific canonical markers (in addition to the cell type-specific canonical markers in Fig 4A) for 12 clusters with corresponding cell types ( Fig 4B). Additionally, a heatmap shows the canonical cell type markers from all the immune cell types along with blood endothelial cells (BECs) and lymphatic endothelial cells (LECs) (S3 Fig).

Detection of Leishmania major transcripts in multiple cell types other than macrophages
We aligned the reads to Leishmania major (LM) reference genome to detect the presence of Leishmania transcripts in 35 different cell type clusters. Interestingly, the differential expression of LM transcripts from scRNA-Seq revealed 20 of 35 cell types have at least one LM transcript associated with that cell. As predicted, we found macrophages are the top immune cell type with about 10% of cells associated with LM transcripts (Fig 5A and 5B). We found at least 2% of other immune cell types are also associated with LM transcripts including resident macrophages, DCs, and neutrophils. At least 1% of the cells in CD4 + Th cells, CD8 + cytotoxic T cells, T regulatory cells, and basophils were also associated with detectable LM transcript ( Fig  5A and 5B). Consistent with previous findings, we found fibroblasts and keratinocytes were also associated with LM transcripts [59]. Surprisingly, we also detected LM transcripts associated with >5% myoblasts, but not myocytes and almost 1% of the BECs. Overall, our data shows that multiple other cell types at the infection site are associated with LM transcripts suggesting other cell types maybe infected with parasites alongside monocytes and macrophages, which are the well-established primary host cells for Leishmania parasites.

scRNA-Seq confirms the immune cell recruitment at the site of L. major infection in vivo
Murine L. major infection leads to the recruitment of immune cells such as neutrophils, inflammatory monocytes, and monocyte-derived macrophages to the site of infection. Particularly, inflammatory monocytes and both tissue-resident macrophages and monocyte-derived macrophages can serve as a replicative niche for the pathogen, as well as play a role in parasite control by killing the pathogen. BECs mediate immune cell recruitment to the infected and inflamed tissue and LECs promote immune cell migration away from the infected skin.
Therefore, we speculate that immune cells and ECs participate in parasite control and/or immunopathology during CL. As a result, the remainder of the study focuses on 7 immune cell types along with the BEC and LEC clusters. Our UMAP projection displays 13,034 cells in naive ears and 13,524 cells in the infected ears ( Fig 6A). Consistent with previous findings in CL, the UMAP plot confirms a significant recruitment of various immune cell types such as inflammatory monocytes, neutrophils, macrophages, dendritic cells, NK cells, and CD4 + and CD8 + T cells in the infected ears that are seen at higher frequencies compared to naive controls ( Fig 6A). Concordant with our scRNA-Seq results, flow cytometric analysis detected a significant increase in the frequency and cell number of macrophages ( Fig 6B and 6C), Ly6C + inflammatory monocytes (Fig 6D and 6E) and neutrophils (Fig 6F and 6G) in infected ears compared to naive controls, while no significant alterations in the BEC or LEC populations were observed (Fig 6H-6J). Altogether, these data confirm the enhanced immune cell migration during L. major infection and transcriptional changes within these individual cell types were investigated.

Differential gene expression of immune cell types during L. major infection
To explore the transcriptional changes in a cell type-specific manner, the DEGs were compared between infected and naive mice within an individual cell type following L. major infection. A volcano plot showing DEGs for macrophages, resident macrophages, inflammatory monocytes, and neutrophils reveals several markers indicative of an increase in myeloid cells in leishmanial lesions (Fig 7A-7D). For instance, transcripts commonly elevated within the top 10 DEGs in myeloid cells and dendritic cells (DCs) include B2m, H2-K1, Gbp2, ligp1, whereas multiple ribosomal proteins showed reduced transcript abundance within the top 10 DEGs among myeloid cells and DCs (Tables 2 and 3 (Table 2A), resident macrophages (Gbp4, Gbp8, Gbp2) (Table 2B), inflammatory monocytes (Gbp2, Gbp5, Gbp7, Gbp3) (Table 3A), and DCs (Gbp2) (S4 Fig). Many of the transcriptomic differences detected in myeloid cells, such as elevated GBPs, were also found in EC populations. For instance, BECs expressed increased Gbp4 and Gbp2, and LECs expressed increased Gbp4, Gbp2, and Gbp7 upon infection   (Table 4A- In addition, many chemokines were differentially modulated in myeloid cells and ECs with L. major infection (Fig 7). Specifically, Cxcl9 was significantly elevated in macrophages, resident macrophages, inflammatory monocytes, DCs, BECs, and LECs following L. major infection (Figs 7 and S4). In contrast, we found significant downregulation of Ccl24 in resident macrophages, Cxcl12 in BECs, and Ccl17 in DCs following L. major infection (Figs 7B and 7E and S4). Also important in immune cell recruitment, selectins (Selp, Sele) and adhesion molecules (Vcam1) were significantly upregulated in BECs with infection, while tight junction molecules like Cldn5 were downregulated (Fig 7E). Of note, known canonical markers were significantly elevated with L. major infection including Arg1, Nos2, and Pla2g7 in macrophages, Fcgr4, C3, and Ccl8 in resident macrophages, and Ifitm1, Syngr2, Cd200, Ccr7 in DCs; the complete list of DEGs that are enriched during L. major infection from these and other cell type clusters such as fibroblasts, keratinocytes, chondrocytes, sebaceous glands, basophils, upper hair follicle cells, pericytes, schwann cells, mast cells, myocytes and myoblasts can be found in Gene Expression Omnibus database with the GEO accession number-GSE181720.

Characterization of upstream gene regulators and canonical pathways during L. major infection in vivo
Next, IPA analysis was performed to define the transcriptomic signature for each individual cell type at the site of L. major infection. IPA analysis of our scRNA-Seq data revealed several known and unknown canonical pathways, upstream regulators, and disease-based functional networks. Here, we present the transcripts that are significantly altered (adj. p value < 0.05) in macrophages, BECs, and LECs from infected ears compared to naive controls.

Canonical pathways
IPA analysis highlighted an unknown role for eukaryotic translation initiation factor 2 (EIF2) signaling which includes large ribosomal proteins (Rpl) and small ribosomal proteins (Rps) that had a significantly lower RNA abundance with infection compared to naive animals ( Fig  8). Remarkably, the EIF2 pathway was the top downregulated pathway amongst multiple cell types including macrophages (Ranked as 1 amongst 377), BECs (Ranked as 1 amongst 257), and LECs (Ranked as 1 amongst 145) (Fig 8A and 8C and 8E). This was followed by an involvement of mTOR signalling and eIF4/p70S6K signalling in macrophages, BECs, and LECs from infected ears. Alongside the IPA pathway results, the expression of individual transcripts for each of the corresponding pathways is provided for macrophages (Fig 8B), BECs (Fig 8D), and LECs ( Fig 8F). These data show the top 10 transcripts in the EIF2 signalling PLOS NEGLECTED TROPICAL DISEASES pathway in macrophages, BECs, and LECs, which includes many subunits of ribosomal proteins and the RNA abundance of these transcripts is mostly decreased in infected ears compared to naive controls (Fig 8B and 8D and 8F). In addition, we validated the relative expression of some of the Rpl and Rps transcripts by RT-PCR and found the significantly lowered expression of Rpl4, Rpl12, Rps18, Rps19 during the infection compared to naïve animals (Fig 9A-9D). In contrast, the IPA analysis revealed the antigen presentation pathway was increased with infection, and this pathway was also a common feature of infection being the top elevated pathway for macrophages, BECs, and LECs (Fig 8B and 8D and 8F). The antigen presentation pathway was enriched in transcripts such as B2m, Cd74, H2-K1, H2-Aa, H2-Ab1, H2-Eb1 that were elevated with L. major infection in macrophages, BECs, and LECs. In addition to macrophages, BECs, and LECs, we also noted that an involvement of mTOR signalling pathway is consistent among the other immune cell types apart of this study such as inflammatory monocytes, DCs, and CD4 + T cells, as mTOR signalling is top 5 in the list of pathways (S6 Fig). In summary, EIF2 signaling is the top downregulated pathway in BECs, LECs, macrophages, as well as other immune cell types from infected ears.

Discussion
We conducted a comprehensive high-resolution transcriptomic analysis using both bulk and scRNA-Seq approaches to discover the global changes in the transcript expression that occurs following L. major infection in vivo. Through our bulk RNA-Seq analyses, we identified many differentially regulated novel transcripts in immune compartments that are related to host immune response pathways. We found significant enrichment of DEGs in the antigen processing and presentation pathway following L. major infection. Specifically, our data indicate that L. major infection upregulates many MHC molecules belonging to the antigen processing and presentation pathway along with inflammatory cytokines such as IFNγ following L. major infection. These findings are consistent with the well-established role of the Th1 immune response in parasite control; we confirm antigen presentation is a process associated with the host response to Leishmania infection [60][61][62][63][64][65]. Additionally, we noted the enrichment of DEGs specific for other host immune response pathways such as chemokine signalling, cell adhesion molecules, and cytokine-cytokine receptor interactions in infected ears. In contrast, our bulk RNA-Seq results revealed DEGs associated with the ribosomal pathway were downregulated following L. major infection. While the importance of the decrease in RNA abundance of transcripts that encode 40S and 60S ribosomal subunits has not been studied in Leishmania, these findings suggest cells at the site of infection are actively controlling translation and/or ribogenesis or undergoing a stress response similar to bacterial, viral infection and cancer [66][67][68]. Altogether, our bulk RNA-Seq results confirm a known role for the importance of antigen presentation and highlight an unknown feature of downregulating ribosomal subunits in CL.

Host transcriptomic changes following Leishmania major infection
Our scRNA-Seq analysis revealed many novel transcripts from various cell types that remain largely unexplored at the site of L. major infection. Through our scRNA-Seq data, we defined 35 distinct cell type populations by using canonical markers specific for different cell types including both resident and recruited cells following L. major infection. In agreement with previously published results [11,69], our scRNA-Seq analysis confirmed a significant increase in various immune cell types at the site of infection with recruited myeloid cells such as neutrophils, inflammatory monocytes, and monocyte-derived macrophages in lesions ( Fig  6). The transcriptional signature of interferon-induced GTPases like GBPs was significantly upregulated in macrophages, resident macrophages, inflammatory monocytes, DCs, BECs, and LECs following L. major infection. These data suggest GBPs may play a protective role in both immune and nonimmune cells during L. major infection. GBPs are involved in controlling intracellular pathogen replication and specifically mediating the protection against intracellular pathogens such as Listeria monocytogenes and Mycobacterium bovis [70] [71]. Importantly, mice deficient in GBP genes are more susceptible to Toxoplasma gondii [72]. Moreover, GBPs restrict L. donovani growth in nonphagocytic cells such as murine embryonic fibroblasts in an IFNγ-dependent manner [73].
We detected several transcripts for chemokines in myeloid cells that are elevated with infection that may mediate myeloid cell accumulation at the site of infection. Similar to the transcriptomic findings of Carneiro et al., we detected increased RNA abundance of Ccl5, Cxcl9 and Cxcl10 with infection, whereas Ccl2 and Ccl7 was not detected in inflammatory monocytes and macrophage population in our data set [16]. However, unlike Carneiro et al., we did not separate infected from bystander inflammatory monocytes in our studies. CXCL9 and CXCL10 are monocyte-recruiting chemokines produced in response to IFNγ further supporting a potential role for Th1 immunity in recruiting permissive Ly6C + CCR2 + monocytes during Leishmania infection [16]. Our results revealed a significant downregulation of Ccl24 transcript in resident macrophages from infected ears compared to naive controls. A recent study demonstrated dermal tissue-resident macrophages shift toward a pro-inflammatory state in L. major-infected mice lacking IL-4/IL-13 from eosinophils, which is mainly regulated by CCL24 from resident macrophages [74]. Also, potentially mediating myeloid cell migration, we found BECs in infected skin express elevated transcripts for selectins (Sele, Selp) and adhesion molecules (Vcam1), while concomitantly downregulating transcripts responsible for junctional stability. Mice doubly deficient in E-and P-selectin develop significantly less inflammation following L. major infection [75]. Taken with our findings, these data suggest BECs play an active role during CL to recruit immune cells into the site of infection. Indeed, these scRNA-Seq data also correlate with our bulk RNA-Seq results which indicate a role of leukocyte trans-endothelial migration pathway in infected ears when compared to naive controls. Overall, our scRNA-Seq data reveals cells in leishmanial lesions exist in pro-inflammatory environment, but the actual host protective role of individual DEGs within each cell type requires further investigation.
We detected the presence of LM transcripts associated with multiple cell types including myeloid cells like macrophages, inflammatory monocytes, neutrophils, and DCs which are all known to harbor parasites [59,76,77]. We also detected LM transcripts associated with stromal cells such as fibroblasts and keratinocytes, which can be infected by parasites [78,79]. Surprisingly, we also found evidence of parasite transcripts associated with myoblasts, chondrocytes, and BECs which has not previously been reported. Importantly, the molecular techniques used in these studies cannot differentiate between living and dead parasites and we can only report cells associated with parasite transcripts. Therefore, subsequent studies will need to determine if parasites in other non-myeloid cells like myoblasts and BECs are viable and capable of replication and infection. Others have hypothesized Leishmania parasites might evade the host immune responses by seeking shelter in different non-macrophage cell types including fibroblasts and keratinocytes in addition to infection in myeloid cells which exhibit a more robust pro-inflammatory response against Leishmania parasites [80]. However, the presence of L. major transcripts from non-myeloid stromal cell types like myoblasts and BECs needs to be further explored to determine whether these cells can serve as safe havens for L. major parasites during chronic infection or provide a conduit for metastasis. It should be noted there are limitations to these findings as the current analysis cannot distinguish between cells with internalized parasites compared to cells with parasites that are co-purified with cells given free parasites are present in tissues and parasites could be attached to the outside surface of cells. Regardless, these results shed light on cell types previously not thought to harbor Leishmania parasites in the skin and provide an opportunity for new investigation. Lastly, future studies will compare the transcriptomic profiles of cells that are associated with LM transcripts with cells that are not are associated with LM transcripts.
Inflammatory monocytes and macrophages serve as replicative niches for the parasite as well as the major cell types responsible for parasite control. BECs play a crucial role in regulating immune cell entry into inflamed tissues and LECs participate in immune cell migration out of the lesions. As a result of the critical roles of these cell types during L. major infection, we focused on identifying the upstream regulators and canonical pathways for macrophages, BECs, and LECs. Surprisingly, the antigen processing and presentation pathway and EIF2 signaling were the most significant pathways in all three cell types (macrophages, BECs, and LECs) within infected ears. The antigen processing and presentation showed a positive activation z score indicating overall upregulation of the pathway, which is consistent with findings in human CL lesions [81]. Immunoproteasomes play a critical role in the immune response by degrading intracellular proteins to generate MHCI epitopes for effective antigen presentation [82]. While the increased expression of immunoproteasome genes in human lesions caused by L. braziliensis infection has been reported [81], our results reveal for the first time that LECs express higher levels of transcripts for immunoproteasomes (psmb8 and psmb9), as well as transcripts involved in the antigen presentation pathway (Tap1 and Tapbp) following L. major infection. These data suggest that ECs, and specifically LECs, may play an unknown role in antigen presentation during L. major infection, similar to viral infection and vaccination [82,83].
The EIF2 signaling pathway had the highest negative activation z score of all the pathways indicating overall downregulation during CL. To our knowledge, this is the first study highlighting EIF2 signaling as a novel candidate pathway for leishmaniasis. Eukaryotic initiation factor-2 (EIF2) is a GTP-binding protein, which initiates protein translation by delivering charged initiator met-tRNA onto the ribosome [84]. Upon subject to infection-induced cellular stress, EIF2 plays a significant role in attenuating translation initiation by phosphorylation of the alpha subunit of eIF2 leading to immediate shut-off of translation and activation of stress response genes [84]. Phosphorylation of eIF2α plays as a rate limiting step as it reduces active eIF2-GTP levels at translation initiation which ultimately results in a global reduction of protein synthesis [84,85]. Our data indicate that impaired EIF2 signaling is linked to the downregulation of many ribosomal subunits in macrophages, BECs, and LEC in L. major-infected ears (Fig 8). The known phenomenon of "protein shut off" has been well described at the molecular level for some viruses [86], but there is little to no evidence documenting this phenomenon during L. major infection. Generally, the enhanced host response to viral and bacterial infections depends on the upregulation of EIF2-mediated translational control, thereby reducing general protein synthesis [67,68,86]. However, we found EIF2 signaling was downregulated with L. major infection. Previously, it was demonstrated that L. major promotes its survival by downregulating macrophage protein synthesis, which is mainly mediated by host translation repressor 4E-BP1 [87]. Given the robust production of cytokines and chemokines in leishmanial lesions, global protein synthesis does not seem impacted by L. major infection in vivo, but careful analysis of EIF2 signaling has not been performed to date.
Alongside EIF2 signalling, the involvement of mTOR and eIF4/p70S6K signalling in macrophages, BECs, LECs, DCs and inflammatory monocytes in infected mice hints at an important role for hypoxia-induced oxidative stress at the site of L. major infected ears. Our results indicate the activation of metabolic gene targets like hypoxia-inducible factor (HIF-1α) in macrophages, which is consistent with our previous results and work by the Jantsch lab showing leishmanial lesions are hypoxic [9,[88][89][90]. However, the consequence of hypoxia-induced oxidative stress and/or HIF-1α signaling is complex, and appears to be context dependent based on the cell type or Leishmania species as HIF-1α hampers DC function during VL, but enhances macrophage effector responses in both parasite control and orchestrating vascular remodeling in CL [9,10,[89][90][91][92]. We speculate EIF2 signaling is downregulated as part of the stress response, potentially from hypoxia, but mechanism by which EIF2 signaling is impaired in multiple cell types and how that contributes to pathogen control or the pathogenesis of disease in CL is unknown. A caveat to our studies is that we only examined a single time point at 4 weeks p.i. following L. major inoculation. Over the course of infection, parasite burdens and lesion severity changes with time. The lesion microenvironment is dynamic and composed of different cell types at different time points following infection, especially as the inflammatory reaction subsides and the wound healing process takes over. However, these findings will inform subsequent studies examining the transcriptomic profiles of individual cell types at other time points following inoculation such as 1-2 weeks p.i. when the Th1 response is active but lesions are not present, at 8-10 week p.i. when lesions are healing, and at 12 weeks p.i. when lesions have resolved. Immune cells are recruited into lesions but our analysis also does not take into consideration how long a cell has resided in the lesion, being exposed to soluble mediators and the hypoxic microenvironment which could dramatically reprogram the transcriptomic signature of that cell.
Collectively, our transcriptome analysis not only provides the first comprehensive list of gene expression at single-cell resolution, but also highlights a previously unknown role of the highly conserved EIF2 signaling pathway in leishmaniasis. At the same time, our study has some limitations. This work was only performed in a resistant C57BL/6 mouse model, and it is not clear if the transcriptomic findings such as the decrease in EIF2 signaling and ribosomal transcripts are a general feature of CL. Future studies will need to determine if our findings will be replicated in a susceptible BALB/c mouse infected with L. major. Similarly, this study was only performed with L. major and it remains unknown if other healing and nonhealing CL species like L. mexicana or L. amazonensis will yield similar results. Regardless, future analysis by us and others utilizing these datasets will expand our knowledge on the complex immune and stromal cell networks and pathways participating in the host response to Leishmania infection.

Ethics statement
All animal procedures were performed in accordance with the guidelines of the UAMS Institutional Animal Care and Use Committee (IACUC) under the animal protocol 4013.

Animals
Female C57BL/6NCr mice were purchased from the National Cancer Institute. Mice were housed in the Division of Laboratory Animal Medicine at University of Arkansas for Medical Sciences (UAMS) under pathogen-free conditions and used for experiments between 6 and 8 weeks of age.

Parasite infection in vivo
Leishmania major (WHO/MHOM/IL/80/Friedlin) strain was used. Parasites were maintained in vitro in Schneider's Drosophila medium (Gibco) supplemented with 20% heat-inactivated FBS (Invitrogen), 2 mM L-glutamine (Sigma), 100 U/mL penicillin, and 100 mg/mL streptomycin (Sigma). Metacyclic stationary phase promastigotes were isolated from 4-5 day cultures by Ficoll density gradient separation (Sigma), as previously described [93]. For ear dermal infections, 2×10 6 L. major promastigote parasites in 10 μL PBS (Gibco) were injected intradermally into the ear. Lesion development was monitored weekly by measuring ear thickness and lesion area with a caliper. Lesion volume was calculated. At 4-week post infection, ears were excised and enzymatically digested using 0.25 mg/mL Liberase (Roche) and 10 mg/mL DNase I (Sigma) in incomplete RPMI 1640 (Gibco) for 90 min at 37˚C. After digesting, ears were minced manually to obtain cellular content in single cell suspensions. Parasite burdens were determined by limiting dilution assays (LDAs), as previously described [94].

Bulk RNA-Seq: Sample preparation
The total RNA was isolated from the cell lysate of naive and infected ears by using Qiagen's RNeasy plus mini kit according to the manufacturer's instructions. The CTPR Genomics and Bioinformatics Core at the Arkansas Children's Research Institute (ACRI) prepared sequencing libraries from RNA samples by use of the Illumina TruSeq Stranded mRNA Sample Preparation Kit v2. for sequencing on the NextSeq 500 platform using Illumina reagents. The quality and quantity of input RNA was determined using the Advanced Analytical Fragment Analyzer (AATI) and Qubit (Life Technologies) instruments, respectively. All samples with RQN (RNA quality number) values of 8.0 or above were processed for sequencing. Sequencing libraries were prepared by use of the TruSeq Stranded mRNA Sample Prep Kit (Illumina). Briefly, total RNA (500 ng) was subjected to polyA selection, then chemically fragmented and converted to single-stranded cDNA using random hexamer primed reverse transcription. The second strand was generated to create double-stranded cDNA, followed by fragment end repair and addition of a single A base on each end of the cDNA. Adapters were ligated to each fragment end to enable attachment to the sequencing flow cell. The adapters also contain unique index sequences that allow the libraries from different samples to be pooled and individually identified during downstream analysis. Library DNA was PCR amplified and enriched for fragments containing adapters at each end to create the final cDNA sequencing library. Libraries were validated on the Fragment Analyzer for fragment size and quantified by use of a Qubit fluorometer. Equal amounts of each library was pooled for sequencing on the NextSeq 500 platform using a high output flow cell to generate approximately 25 million 75 base reads per sample.
Genes with unique Entrez IDs and a minimum of~2 counts-per-million (CPM) in 4 or more samples were selected for statistical testing. This was followed by scaling normalization using the trimmed mean of M-values (TMM) method [98] to correct for compositional differences between sample libraries. Differential expression between naive and infected ears was evaluated using limma voomWithQualityWeights [99] with empirical bayes smoothing. Genes with Benjamini & Hochberg [100] adjusted p-values � 0.05 and absolute fold-changes � 1.5 were considered significant.
Gene Set Enrichment Analysis (GSEA) was carried out using Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway databases and for each KEGG pathway, a p-value was calculated using hypergeometric test. Cut-off of both p < 0.05 and adjusted p-value/FDR value < 0.05 was applied to identify enriched KEGG pathways. DEGs that are more than 1.5-fold in L. major-infected ears relative to uninfected controls were used as input, with upregulated and downregulated genes considered separately. Subsequently, the heat maps were generated using these genes with complex Heatmap. All analyses and visualizations were carried out using the statistical computing environment R version 3.6.3, RStudio version 1.2.5042, and Bioconductor version 3.11. The raw data from our bulk RNA-Seq analysis were deposited in Gene Expression Omnibus (GEO accession number-GSE185253).

scRNA-Seq sample preparation
The Genomics and Bioinformatics Core at the Arkansas Children's Research Institute (ACRI) prepared NGS libraries from fresh single-cell suspensions using the 10X Genomics NextGEM 3' assay for sequencing on the NextSeq 500 platform using Illumina SBS reagents. The quantity and viability of cells input to the assay was determined using Trypan Blue exclusion under 10X magnification. Library quality was assessed with the Advanced Analytical Fragment Analyzer (Agilent) and Qubit (Life Technologies) instruments.

scRNA-Seq data analysis
Demultiplexed fastq files generated by the UAMS Genomics Core were analyzed with the 10X Genomics Cell Ranger alignment and gene counting software, a self-contained scRNA-Seq pipeline developed by 10X Genomics. The reads are aligned to the mm10 and Leishmania major reference transcriptomes using STAR and transcript counts are generated [96,101].
The raw counts generated by cellranger count were further processed using the R package Seurat [102,103]. Low quality cells, potential cell doublets, and cells with high percentage of mitochondrial genes were filtered out of the data. We filtered cells that have unique feature counts over more the 75 th percentile plus 1.5 times the interquartile range (IQR) or less than the 25 th percentile minus 1.5 time the IQR. Additionally, we filtered cells with mitochondrial counts falling outside the same range with respect to mitochondrial gene percentage. Following filtering the counts all 8 sequencing runs were merged. The counts are then normalized using the LogNormalize method, which normalizes the feature expression measurements for each cell by the total expression, multiplies this by a scale factor (10,000 by default), and logtransforms the result. Next, the 2000 highest variable features are selected. The data is then scaled, and linear regression is included to remove variation associated with percent mitochondria and cell cycle status. Principle component analysis (PCA) is performed on the scaled data. A JackStraw procedure was implemented to determine the significant PCA components that have a strong enrichment of low p-value features.
A graph-based clustering approach is applied [104]. Briefly, these methods embed cells in a graph structure-for example a K-nearest neighbor (KNN) graph, with edges drawn between cells with similar feature expression patterns, and then attempt to partition this graph into highly interconnected 'quasi-cliques' or 'communities'. t-distributed stochastic neighbor embedding (tSNE) and Uniform Manifold Approximation and Projection (UMAP) [105] are non-linear dimensional reduction techniques used to visualize and explore the results and are performed using Seurat. Seurat FindNeighbors and FindClusters functions are optimized to label clusters based on the visual clustering in the projections. Seurat FindAllMarkers function finds markers that define clusters by differential expression. It identifies positive markers of a single cluster compared to all other cells and outputs the differential expression results. These markers are compared to known markers of expected cell types and results from previous single-cell transcriptome studies in order to assign appropriate cell type labels. Cell type determinations were made by manually expecting these results and some clusters were combined if their expression was found to be similar. Differential expression analysis is performed using MAST, a GLM-framework that treats cellular detection rate as a covariate [106]. The raw data from our scRNA-Seq analysis were deposited in Gene Expression Omnibus (GEO accession number-GSE181720).

Ingenuity Pathway Analysis (IPA)
QIAGEN's Ingenuity Pathway Analysis (IPA, QIAGEN Redwood City, www.qiagen.com/ ingenuity) was utilized to investigate the DEGs at the level of biochemical pathways and molecular functions. We submitted our DEGs to functional analysis with IPA, and IPA provided canonical pathways, diseases and function, and upstream regulators based on the experimentally observed cause-effect relationships related to transcription, expression, activation, molecular modification, etc. Z-score analyses are used to assess the match between observed and predicted up and down regulation patterns allowing for Bayesian scoring of the results.