Transcriptional profiling of macrophages reveals distinct parasite stage-driven signatures during early infection by Leishmania donovani

Macrophages undergo swift changes in mRNA abundance upon pathogen invasion. Herein we describe early remodelling of the macrophage transcriptome during infection by amastigotes or promastigotes of Leishmania donovani. Approximately 10–16% of host mRNAs were differentially modulated in L. donovani-infected macrophages when compared to uninfected controls. This response was partially stage-specific as a third of changes in mRNA abundance were either exclusively driven by one of the parasite forms or significantly different between them. Gene ontology analyses identified categories associated with immune functions (e.g. antigen presentation and leukocyte activation) among significantly downregulated mRNAs during amastigote infection while cytoprotective-related categories (e.g. DNA repair and apoptosis inhibition) were enriched in upregulated transcripts. Interestingly a combination of upregulated (e.g. cellular response to IFNβ) and repressed (e.g. leukocyte activation, chemotaxis) immune-related transcripts were overrepresented in the promastigote-infected dataset. In addition, Ingenuity Pathway Analysis (IPA) associated specific mRNA subsets with a number of upstream transcriptional regulators predicted to be modulated in macrophages infected with L. donovani amastigotes (e.g. STAT1 inhibition) or promastigotes (e.g. NRF2, IRF3, and IRF7 activation). Overall, our results indicate that early parasite stage-driven transcriptional remodelling in macrophages contributes to orchestrate both protective and deleterious host cell responses during L. donovani infection.


Results
Infection with L. donovani amastigotes or promastigotes promotes early changes in the mRNA pool of the host cell. To compare the early effects of the two life stages L. donovani in the mature mRNA pool of the host cell, total cytosolic mRNA extracts from bone marrow-derived macrophage (BMDM) cultures infected with amastigote (AMA) or promastigote (PRO) parasites for 6 h were subjected to RNAseq and compared to non-infected controls (CTR) (Fig. 1A). As shown by a projection of a principal component analysis, infection appears to be the main source of variation (PC1 = 37.4%) between the different datasets followed by  Table S1). These data indicate that infection by either amastigotes or promastigotes of L. donovani leads to early reprograming of the mRNA content of the host cell.
Early transcriptional changes in macrophages infected with L. donovani amastigotes are associated with the inhibition of cell death and immune functions. Gene Ontology (GO) hierarchical clustering analysis was carried out to determine whether subsets of mRNAs encoding functionally related proteins are selectively modulated in BMDMs upon infection with L. donovani amastigotes ( Fig. 2A and Table S2). Enrichment of functional categories related to regulation of gene expression, positive regulation of DNA repair, and negative regulation of apoptosis and protein modification was detected in the AMA-upregulated dataset ( Fig. 2A Fig. 3B and Table S2). In parallel, mRNAs encoding proteins associated with cell survival (Hmox1, Hsp90ab1, Optn, Wfs1), iron transport (Slc11a2, Slc25a37, Slc39a14, Slc40a1) and redox homeostasis (Cat, Gclm, Gsr, Prdx1, Sod2, Txnrd1) were also overrepresented in the upregulated dataset ( Fig. 3A upper panel, Fig. 3B and Table S2). In line with previous observations 6,24 , an increase in the abundance of a group of transcripts associated with lipid metabolism was detected in the PRO data set (Cd36, Lrp12, Lpl, Acsl1, Fabp4) (Table S1). In contrast, GO categories related to cell death (Casp2, Casp6, Cradd, Dfna5, Dusp6, Mef2c, Rassf2, Sarm1) and immune functions such as leukocyte activation (Clec4a2, Dock8, Gpr183,  Table S2). These data indicate that L. donovani promastigote infection elicits a transcriptome-wide response in macrophages that results in the upregulation of lipid metabolism, the concomitant expression of activating and inhibitory immune mediators, and the inhibition of cell death and antigen presentation.
Parasite stage-specific modulation of the host cell transcriptome during L. donovani infection. Anota2seq Table S2). The PRO enhanced UP subset showed an overrepresentation of apoptosis inhibitors (Bcl2a1d, Gbe1, Gclc, Hmox1, Il1rn, Plk2, Serpinb9) (Fig. 5A, right panel, and Table S2). Consistent with this, the activation of a transcriptional regulatory network leading to the inhibition of cell death was identified by Ingenuity Pathway Analysis (IPA) in the PRO upregulated subset (Fig. S1). Unlike the PRO enhanced transcripts, no GO categories were enriched in the AMA enhanced subset   (Table S3). Some upstream regulators were common between both parasite stages (MYC, KLF4, and SMAD3) albeit with variations in the number and/or identity of downstream targets in each type of infection ( Fig. 6A left panel and Table S3). Others were predicted to be activated only by amastigotes (YY1, WDR5, and TP73) or promastigotes (NFE2L2, IRF7, IRF3, EPAS1, SPI1, NFATC2, ATF4, IFI16, CEBPB, CREB1, SP1, FOXO1, and FOS) (Fig. 6A left panel and Table S3). As expected, transcriptional regulators predicted to be activated upon L. donovani infection showed high percentages of associated upregulated mRNAs (Fig. 6B). In agreement with predicted induction of NFE2L2 (a.k.a. NRF2)-dependent transcriptional programs in BMDMs infected with L. donovani promastigotes (i.e. 63 genes) (Fig. 6B right panel, Fig. S2A and Table S3), NRF2-mediated Oxidative Stress Response was identified by IPA as one of the top networks to be activated by the promastigote stage (Fig. S2B). In addition, a small group  Table S3). Of note, SIRT1 was predicted to be activated in the amastigote-infected dataset (Fig. 6A left panel and Table S3) whereas the opposite was observed during infection with the promastigote stage (Fig. 6A right panel and Table S3), as previously reported 25 . These data hint at the involvement of a complex regulatory network affecting the abundance of functional subsets of mRNAs in BMDMs infected with L. donovani amastigotes or promastigotes.

Discussion
Early remodelling of the macrophage transcriptome has been reported to be pathogen-specific during bacterial and parasitic infections 17,24,[26][27][28] . Transcriptome-wide analyses of macrophages infected with L. donovani have mainly been described at ≥ 12 h post-infection 18,20,21,[29][30][31] , thereby omitting an earlier timeframe during which numerous molecular and cellular changes occurring within infected macrophages 8,9,11,32 could trigger, or be elicited by, selective reprogramming of the host transcriptome. Herein, using RNAseq, we describe rapid changes in the levels of mRNAs of primary murine macrophages infected with L. donovani amastigotes and promastigotes. Distinct transcriptional signatures were identified in macrophages infected with each parasite stage. A marked inhibition of mRNAs encoding proteins related to different immune functions was found in the amastigote-infected dataset whereas a combination of activating and inhibitory immune modulators was observed in promastigote-infected macrophages. Additionally, our in silico analyses identified host mRNA signatures in the up-and downregulated datasets that appear to be under the control of parasite-stage driven networks of transcription factors. These observations indicate that amastigotes and promastigotes of L. donovani elicit a complex transcriptome-wide reprogramming in infected macrophages that includes both parasite stage-specific and commonly regulated mRNA subsets. www.nature.com/scientificreports/ Leishmania donovani amastigote-driven changes in macrophage gene expression have been documented at ≥ 24 h post-infection 18,30,33 . Herein, we provide evidence that L. donovani amastigote infection leads to a vast remodelling of the macrophage transcriptome as early as 6 h post-infection. Among the downregulated targets, we found an enrichment in mRNAs encoding proteins related to several macrophage immune functions. IPA predicted that some of these changes are dependent on the inhibition of transcription factor STAT1. In this regard, Matte and Descoteaux previously reported that L. donovani amastigotes prevent STAT1 nuclear import and pro-inflammatory gene expression (i.e., Nos2 and Irf1) in BMDMs stimulated with IFNγ 13 . In addition, a transcriptomic study carried out in splenic macrophages revealed that these cells become insensitive to IFNγ during experimental VL despite a strong pro-inflammatory environment in the spleen 30 . Hence, it is plausible that early blockade of STAT1-dependent transcriptional programs in macrophages infected by L. donovani amastigotes has a negative effect in IFNγ-mediated microbicidal and immune host responses at later stages of the disease. Further investigation is required to shed light on this matter.
Infection of macrophages results in an oxidative burst response that involves the production of potent microbicidal effectors such as reactive oxygen and nitrogen species 34 . However, the antimicrobial oxidative stress response can also compromise macrophage DNA integrity and lead to the activation of apoptotic signals 35 40 . Furthermore, upregulation of NRF2 activity contributed to promote parasite persistence during L. guyanensis infection by limiting inflammation 40 . In addition, NRF2-dependent increase in heme oxygenase 1 (HO-1) and ATF3 upon L. donovani infection was critical in dampening macrophage oxidative burst and proinflammatory cytokine expression as part of a parasite survival strategy 15 . Thus, our data showing an enrichment of transcripts associated with the activation of an NRF2dependent antioxidant response in promastigote-infected BMDMs suggest that targeting this regulatory node could be a therapeutic approach to combat VL.
Mounting evidence indicates that specific and abundant changes in the transcriptional landscape of macrophages occur with 1-4 h post-infection with promastigotes of different Leishmania species (L. major, L. amazonensis, L. chagasi) [4][5][6][7] . For example, microarray data from BMDMs infected with L. infantum (syn. chagasi) promastigotes for 4 h revealed a marked inhibition of inflammatory transcripts that was concomitant with the upregulation of multiple anti-inflammatory mediators such as TGF-β 7 , a disease severity marker during VL 2 . Even though we did not identify Tgfb1 in the subset of transcripts upregulated in response to early infection with L. donovani promastigotes, we recently described eIF4A-dependent increase in Tgfb1 mRNA translation efficiency in BMDMs infected with L. donovani promastigotes and amastigotes for 6 h 41 . Thus, different VLcausing Leishmania spp. (L. infantum and L. donovani) can lead to similar phenotypes in macrophages, such as rapid production of TGF-β, through different regulatory mechanisms of gene expression.
Our IPA and GO analyses identified a transcriptional signature characterized by early induction of pro-and anti-inflammatory genes in macrophages infected with L. donovani promastigotes. These data are in line with previous reports on early reprogramming of the host cell transcriptome by promastigotes of L. major and L. amazonensis, two Leishmania species that cause cutaneous leishmaniasis (CL). A common feature of this type of signature appears to be the upregulation of the pro-inflammatory gene Tnf (Fig. 3) [4][5][6] . TNF levels have been associated with early recruitment of immune cells, including potential host cells, at the site of infection 42 . Thus, it is conceivable that both VL-and CL-causing Leishmania species drive rapid Tnf transcription and TNF production by macrophages to favor their own replication.
Global-scale profiling of macrophages identified a transcriptional signature associated with the modulation of lipid metabolism during early infection with L. major promastigotes 6 . This was further characterized by showing cholesterol accumulation and the dynamics of lipid droplet formation in infected macrophages 24 . Our in silico analyses identified a subset of lipid metabolism-related mRNAs upregulated in the L. donovani promastigoteinfected data set. Consistent with this, alterations in lipid metabolism have been reported in patients diagnosed with VL 43 . Hence, our data along previous studies indicate that early transcriptional changes triggered by CLand VL-causing Leishmania species contribute to reprogramming lipid metabolism of infected macrophages.
Recently, a transcriptomic analysis of macrophages infected with L. donovani promastigotes identified HIF-1α as a negative regulator of the parasite-promoting BNIP3/mTOR/SREBP-1c lipogenesis axis 23 . In parallel, the induction of a transcriptional signature associated with glutamine metabolism was found to be pivotal in VL pathogenesis with a therapeutic potential in synergy with miltefosine treatment 22 . Both studies performed RNAseq on macrophages infected with L. donovani promastigotes for 6 h and, although identified transcripts were validated in vivo and in vitro, the global transcriptional response of infected macrophages compared to uninfected controls was not analyzed 22,23 . Even though we did not find an enrichment of HIF-1αdependent transcripts in our dataset, we detected an increase in Bnip3, a transcriptional target of HIF-1α, as previously reported 23 . Similarly, our IPA and GO analyses did not find an enrichment of transcripts associated www.nature.com/scientificreports/ with glutamine metabolism; however, mRNAs encoding subunits of glutamate-cysteine ligase, a key enzyme in glutathione synthesis and glutamine usage 44 , were upregulated in infected datasets when compared to uninfected controls (i.e. Gclm in PRO upregulated, and Gclc in PRO and AMA upregulated). In sum, data generated by others and by us indicate that regulation of host cell metabolism is at least in part dependent on parasite-driven transcriptional changes induced by both life stages of L. donovani early during infection.
In line with subversion of macrophage immune functions by L. donovani promastigotes 1 , we identified a number of mRNAs encoding immune inhibitors in the upregulated promastigote-infected dataset, including Cd274 (a.k.a. PDL1), Socs1, and Cd200. PDL1 and its receptor PD1 constitute an important inhibitory axis for T cell activity, and antibody therapy against PD1 has proven successful against numerous malignancies 45 . Notably, the PD1/PDL1 axis was recently identified to play an important role in vivo during VL and immunotherapy against PD1 was effective in hampering parasite burden and pathogenesis 31 . In addition, early induction of SOCS1, a known antagonist of the proinflammatory JAK1/STAT1 pathway 38,46 , was identified as part of a cellular program to prevent oxidative burst-mediated apoptosis in macrophages infected with L. donovani 47 . Similarly, a swift increase of CD200 in macrophages exposed to L. amazonensis or L. donovani infection was described as a strategy to favor parasite proliferation [48][49][50] . Interestingly, immune blockade of CD200 led to an increase in proinflammatory mediators and parasite elimination capacity of macrophages and T cells, showing its potential as a therapeutic target 49 . Taken together, these reports and our transcriptomic study highlight the early ability of L. donovani promastigotes to limit macrophage antimicrobial responses through the modulation of host mRNA abundance.
IPA identified a transcriptional signature associated with type I interferon responses predicted to be activated via the transcription factors IRF3 and IRF7 in the promastigote-upregulated dataset. By contrast, downregulation of Irf7 mRNA abundance was detected in the transcriptome of amastigote-infected BMDMs. IRF7-dependent parasite elimination was reported in macrophages of the splenic marginal zone during the acute phase of L. donovani amastigote infection in vivo (e.g. 5 to 48 h post-infection) and by a cell line of stromal macrophages in vitro. Although the expression of IRF7 was not modulated in hepatic macrophages during VL, IRF7-defficient mice showed a decreased ability to control parasite burden in the liver 51 . These observations along with transcriptomic data and our in silico analysis suggest that the ability of macrophages to elicit IRF7-dependent antimicrobial transcriptional programs upon L. donovani infection is tissue-and/or parasite-stage specific.
Our group recently described rapid remodeling of the translatome of macrophages infected by promastigotes and amastigotes of L. donovani 41 . Herein, we expanded our findings by analyzing early changes in the abundance of host mRNAs during infection. Comparison of the transcriptome and the translatome of L. donovani-infected BMDMs at 6 h post-infection indicates that in contrast to changes in translation efficiency 41 , modulation of mRNA abundance is, at least in part, parasite stage-specific. It is plausible that differences in lipid composition 52 and protein expression 53 between promastigotes and amastigotes can account for these stage-specific profiles. For example, L. donovani promastigotes exhibit a dense glycocalyx comprised of a variety of potent virulence factors (e.g. lipophosphoglycan (LPG), the protease GP63, etc.) that are mostly absent in amastigotes 11,42 . This in turn can affect the process of parasite internalization due to differential usage of macrophage receptors for phagocytosis 54 leading to distinctive host signaling pathways and transcriptional changes upon infection 1 .
Amastigote-driven changes included the upregulation of transcripts encoding DNA repair modulators while inhibiting those encoding antigen-presenting and macrophage activation factors. Alternatively, promastigoteinfected macrophages showed the upregulation of immune inhibitors as well as an antioxidant transcriptional signature associated to NRF2 activity. However, enrichment of transcripts associated with IRF3 and IRF7 suggests that macrophages activate antimicrobial pathways upon L. donovani promastigote infection. Interestingly, mRNAs encoding proteins associated with DNA damage-sensing or DNA repair, apoptosis inhibition and mRNA metabolism were upregulated via changes in abundance (Figs. 2, 3 and 4) and translation efficiency 41 . A similar dual effect was observed on a number of downregulated immune-related transcripts (e.g. antigen presentation, leukocyte activation, etc.) (Figs. 2, 3 and 4) 41 . In all, previous studies, along with our current findings support the notion that early parasite-driven changes in macrophage gene expression programs are under the control of transcriptional and post-transcriptional regulatory mechanisms that tailor both protective and harmful host cell responses during L. donovani infection.

Materials and methods
Reagents and parasites. Culture media and supplements were purchased from Wisent, Gibco, and Sigma-Aldrich. L. donovani (LV9 strain) amastigotes were isolated from the spleen of infected female Golden Syrian hamsters (Harlan Laboratories) as previously described 13 . L. donovani (LV9 strain) promastigotes were differentiated from freshly isolated amastigotes and were cultured at 26 °C in M199 medium supplemented with FBS (10%), hypoxanthine (100 µM), hemin (5 µM), biopterin (3 µM), biotin (1 µM), penicillin (100 U/mL), and streptomycin (100 μg/mL). Early passage stationary phase promastigotes were used for macrophage infections. Differentiation and infection of bone marrow-derived macrophages. Bone marrow-derived macrophages (BMDMs) were differentiated from bone marrow precursor cells isolated from C57BL/6 mice, as previously described 55 . Briefly, marrow was extracted from bones of the hind legs, red blood cells were lysed, and pro- www.nature.com/scientificreports/ genitor cells were resuspended in BMDM culture medium supplemented with 15% L929 fibroblast-conditioned culture medium (LCCM). Non-adherent cells were collected the following day and were cultured for 7 days in BMDM culture medium supplemented with 30% LCCM with fresh medium replenishment at day 3 of incubation. BMDMs were then collected, viable cells were counted by trypan blue exclusion and plated in 150 mm petri dishes at a density of 2 × 10 5 cells/cm 2 overnight. BMDM cultures were inoculated with L. donovani promastigotes or amastigotes at a multiplicity of infection (MOI) of 10:1 for 6 h, as previously described 56 . Glass coverslips were prepared in parallel and stained with HEMA 3 PROTOCOL to assess the rate of infection according to the manufacturer instructions. Promastigote-and amastigote-infected samples averaged at 92.3% ± 2.5% and 86.8% ± 1.9% of infection respectively. Prior to infection, cells were serum-starved for 2 h.
Cytosolic mRNA extraction. Cytosolic lysates of infected and control BMDMs were prepared for RNA extraction as described 55 . RNA was extracted with QIAzol (QIAGEN) and purified using RNeasy MinElute Cleanup Kit (QIAGEN) according to specifications of the manufacturer. Purity and integrity of RNA was assessed using a Bioanalyzer 2100 with an Eukaryote Total RNA Nano chip (Agilent Technologies).
RNAseq and data processing. RNAseq libraries were generated using the Smart-seq2 method 57 and sequenced by using an Illumina HiSeq2500 instrument with a single-end 51-base sequencing setup from three independent biological replicates for uninfected and L. donovani promastigote-infected BMDMs, and five independent biological replicates for L. donovani amastigote-infected BMDMs. First, RNAseq reads mapping to the reference genome of the Nepalese BPK282A1 strain of L. donovani (txid: 981087) were removed (12.7% and 1.4% mappings on average for promastigotes and amastigotes, respectively). The filtered reads were then mapped to the mouse genome assembly GRCm38 (mm10) using HISAT2 with default settings 58 . Gene expression was quantified using the RPKMforgenes.py script 59 with -fulltranscript -readcount -onlycoding flags from which raw per-gene RNAseq counts were obtained (version last modified 07.02.2014). Genes that had zero counts in all samples were discarded. Annotation of genes was obtained from RefSeq.
RNAseq data analysis using anota2seq. RNAseq counts were normalized within anota2seq using the default TMM-log2 method 60 . Significant changes in mRNA abundance were identified by anota2seq 60 using the default parameters with the following modifications: FDR ≤ 0.05; apvEff > log 2 (2.0). In anota2seq, the number of contrasts per analysis equals n-1 being n the number of conditions (i.e. CTR, Ld AMA, Ld PRO). In analysis one, infections were contrasted to the uninfected control (i.e. Ld PRO versus CTR and Ld AMA versus CTR); in analysis two, cells infected by different parasite stages were compared together and an additional contrast was included to complete the anota2seq parameters (i.e. Ld PRO versus Ld AMA and Ld PRO versus CTR). Identifiers for genes which cannot be distinguished based on their high sequence similarity (also reported by RPKMforgenes.py), were excluded from downstream analyses.
Gene ontology analyses. Gene ontology analyses were performed using the PANTHER tool 61  Ingenuity pathway analysis. Enrichment of transcripts showing differential abundance in specific functional networks was determined using Ingenuity Pathway analysis (IPA; QIAGEN) by comparing anota2seqregulated gene sets against the entire sequenced datasets 62 . Within the IPA application, statistical significance was calculated using a right-tailed Fisher Exact test and p-values were adjusted for multiple hypothesis testing using the Benjamini-Hochberg method to arrive at a FDR.
Quantitative RT-PCR. Purified RNA (500 ng) was reverse transcribed using the LunaScript RT SuperMix Kit (New England Biolabs, cat#E3010L). Quantitative PCR was performed with Luna Universal qPCR Master Mix (New England Biolabs, cat#M3003L). Relative quantification was calculated using the comparative Ct method (ΔΔCt) 63 and relative expression was normalized to mouse β-actin. Experiments were performed in independent biological replicates (n = 3); each sample was analyzed in a technical triplicate, the average of which was plotted against the respective conditions used. Primers were designed using NCBI Primer-BLAST (http:// www. ncbi. nlm. nih. gov/ tools/ primer-blast/) (Table S4).