Tissue and host species-specific transcriptional changes in models of experimental visceral leishmaniasis.

Background: Human visceral leishmaniasis, caused by infection with Leishmania donovani or L. infantum, is a potentially fatal disease affecting 50,000-90,000 people yearly in 75 disease endemic countries, with more than 20,000 deaths reported. Experimental models of infection play a major role in understanding parasite biology, host-pathogen interaction, disease pathogenesis, and parasite transmission. In addition, they have an essential role in the identification and pre-clinical evaluation of new drugs and vaccines. However, our understanding of these models remains fragmentary. Although the immune response to Leishmania donovani infection in mice has been extensively characterized, transcriptomic analysis capturing the tissue-specific evolution of disease has yet to be reported.
Methods: We provide an analysis of the transcriptome of spleen, liver and peripheral blood of BALB/c mice infected with L. donovani. Where possible, we compare our data in murine experimental visceral leishmaniasis with transcriptomic data in the public domain obtained from the study of L. donovani-infected hamsters and patients with human visceral leishmaniasis. Digitised whole slide images showing the histopathology in spleen and liver are made available via a dedicated website, www.leishpathnet.org.
Results: Our analysis confirms marked tissue-specific alterations in the transcriptome of infected mice over time and identifies previously unrecognized parallels and differences between murine, hamster and human responses to infection. We show commonality of interferon-regulated genes whilst confirming a greater activation of type 2 immune pathways in infected hamsters compared to mice. Cytokine genes and genes encoding immune checkpoints were markedly tissue specific and dynamic in their expression, and pathways focused on non-immune cells reflected tissue specific immunopathology. Our data also addresses the value of measuring peripheral blood transcriptomics as a potential window into underlying systemic disease.  Conclusions: Our transcriptomic data, coupled with histopathologic analysis of the tissue response, provide an additional resource to underpin future mechanistic studies and to guide clinical research.

Introduction

Of the many diseases associated with infection by the protozoan parasite Leishmania, visceral leishmaniasis (VL) represents one of the most challenging to understand in terms of its pathophysiology. Unlike cutaneous leishmaniasis, where parasite growth and the pathological consequences of the ensuing immune response are largely confined to the site of transmission, VL is characterized by systemic spread of parasites, multi-organ involvement and a systemic response that is fatal without treatment. Systemic cytokine mediated immunosuppression [1, 2], T cell functional defects [3– 7], alterations to the structural integrity of the lymphoid tissue [8– 10] and systemic coagulopathies [11, 12] have all been proposed as mechanisms that promote parasite survival, render the host susceptible to concomitant secondary infections and/or act as negative prognostic markers. Evidence from animal models suggests that tissue specific control over immunity and immunopathology also exists, that may in part reflect extremes within the clinical spectrum of VL [13, 14]. In spite of the broad impact of infection on a range of key physiological processes, a systems-wide appreciation of the underlying immune, metabolic and physiological changes associated with VL has yet to emerge.

Over the past several years, transcriptomic profiling has been adopted as a key methodology for generating an unbiased view of many disease processes. The application of transcriptomic profiling has provided valuable insights into the pathogenesis of many infectious and non-infectious diseases [15– 26], and on the response to drugs, vaccines and other immunotherapeutic interventions [27– 31]. From a clinical perspective, the application of whole blood transcriptomic profiling (WBTP) opened a new era in clinical monitoring of disease, fuelled notably by the work of Chaussebel and others that illustrated its potential to provide insights into systemic disease processes and serve as a clinical tool [19, 32– 34]. Nevertheless, whether WBTP provides information reflecting a subset of systemic events or is an avatar of that response remains to be fully determined. Indeed, studies formally comparing whole blood with the systemic tissue transcriptome are rare [21].

Transcriptomic profiling in leishmaniasis has been limited to date. In human disease, this has largely focused on cutaneous leishmaniasis [35– 37]. A single study has reported on transcriptional changes in the draining lymph node of Sudanese patients with VL before and after treatment with sodium stibogluconate, identifying a potential role for nuclear factor of activated T cells (NFAT)-regulated immune responses in treatment response [38]. One study has examined whole blood transcriptional responses in healthy controls versus asymptomatic and symptomatic VL patients in Brazil [39]. Three studies have evaluated transcriptional changes in the spleen of hamsters infected with L. donovani [40– 42]. Here, we have applied transcriptional profiling to address key aspects of the tissue specific response to L. donovani infection in mice. We describe key differences between immune response and metabolic events occurring in the principal target organs of liver and spleen, as well as in blood, and where possible through publically available data, provide a formal comparison between the splenic response in mice and hamsters and between murine and human whole blood.

Methods

Ethics statement

Experiments were approved by the Animal Welfare and Ethics Review Bodies of the University of York and the LSHTM and the work was performed under UK Home Office license (PPL 60/4377; PPL 70/6997; PPL 70/8207). Mice were killed by cervical dislocation prior to tissue collection, as described below.

Mice and infections

6–8-week-old female BALB/c mice (Charles River, Margate, UK) weighing 20±1 gm and health screened to FELASA 67M standard and maintained under specific pathogen free conditions in individually ventilated cages were used in this study. Leishmania donovani (LV9) parasites were maintained in B6. Rag1 -/- mice and amastigotes prepared following tissue disruption, cell lysis and differential centrifugation, as described elsewhere [43]. 2×10 7 amastigotes in 150 μl RPMI were injected intravenously (i.v.) via the lateral tail vein and without anaesthetic to initiate infection. After infection, mice were allocated to cages of 5 and provided food and water ad libitum. Experiments reported here are derived from two cohorts of mice (designated CRACK-IT_1 and CRACK-IT_2; 20 mice per cohort, n=5 per infection time point and for time-matched uninfected controls). In CRACK-IT_1, mice (n=5) were killed by cervical dislocation at day 15, 17 and 21 post infection. Control naïve mice were killed on the equivalent of day 14 p.i. for logistical reasons. In CRACK-IT_2, infected mice (n=5 per group) were killed at day 36 and 42, alongside a control naïve group. All animals were killed and processed over an approximate time period of 4–6h beginning in the morning. Tissues were aseptically removed post mortem and stored/processed as detailed below. All downstream tissue analysis was performed blinded to group by investigators not involved in animal handling.

Histology and quantitative morphometry

Spleens and livers were embedded in OCT, snap frozen and stored at -80°, before preparing 7 µm cryosections. All subsequent steps were carried out at room temperature. For gross morphology, sections were fixed for 5 minutes in ice-cold acetone and rehydrated by immersion for 2–3 minutes in 90% ethanol, 70% ethanol and finally ddH 2O. Sections were transferred to Harris’s haematoxylin (Sigma-Aldrich, UK) for 5 minutes, then ‘blued’ in running tap water for 5 minutes. Slides were differentiated using 1% acid alcohol for 5 seconds, washed in tap water for 3minutes and stained with 1% w/v eosin (Sigma-Aldrich, UK) in distilled water for 3 minutes. Slides were washed for 3 minutes in ddH 2O and dehydrated through alcohols, cleared in Xylene and mounted in DPX.

For immunohistochemistry, sections were fixed for 5 minutes in ice cold acetone and rehydrated for 2 x 5 minutes in wash buffer (0.05% w/v BSA in PBS) and blocked for 30 minutes in dilution buffer (5% v/v goat serum in wash buffer). After washing, sections were stained for 45 min with AlexaFlour 647-F4/80 or isotype control (BioLegend, USA; 1:500), counterstained with DAPI and mounted in Prolong Gold. Sections were examined using a Zeiss Axio Scan Z1 with a x 20 objective. Granuloma size and distribution was quantified using image analysis software TissueQuest 4.0 (TissueGnostics, Vienna, Austria), blind to treatment group. DAPI was used as the master channel to identify all events within the regions of interest, and clustering of F4/80 + cells was used to identify granulomas. All granulomas within a 1mm 2 region of interest were identified as individual sub regions of interest, and size and cell density was recorded.

Tissue transcriptomics

RNA was isolated from tissue/blood samples and amplified via Agilent low-input Quick Amp labelling kit (Agilent Technologies). Amplified RNA was then assayed with Agilent SurePrint G3 mouse GE 8x60k microarray chips that were scanned with an Agilent C Scanner with Surescan High Resolution Technology (Agilent Technologies). The data were normalized using the percentile shift method to the 75 th percentile. Identification of differentially expressed (DE) genes between infected and naïve samples was performed using the Benjamini and Hochberg false discovery rate (FDR) correction for pairwise comparisons of infected vs. naïve samples [44]. Where multiple probes for a single gene were present, the average expression level was used. This analysis was performed with GeneSpring software (version 9; Agilent) as a standard 5% FDR, with the variances assessed by the software for each t test performed. A 2-fold expression criteria (Log2 FC = 1)was then applied to each gene list. For DE analysis in CRACK-IT_1, a single group of naïve mice were killed at d14 for comparison with samples of infected mice taken at d15, d17 and d21 p.i. For CRACK-IT_2, DE analysis was performed using matched naïve samples for each tissue taken at each time point (see Results).

Gene ontology analysis was performed using GeneSpring (Agilent), and pathway analysis was conducted using Ingenuity Pathway Analysis software (Ingenuity Systems). For the identification of upstream regulators in IPA, we used the z score, which reflects the extent to which known target genes are regulated in the expected direction, whereby positive z scores predict activation of the regulator whereas negative z scores predict inhibition of the regulator. z scores of >2 and <-2 are considered significant. Gene set enrichment analysis (GSEA) [34] was performed to identify enriched gene sets associated with each phenotype (i.e. infected at each time point vs. naïve). Data were collapsed to gene set against an Agilent mouse array chip file (~24K genes) and parameters set for GSEA were: permutations =1000; permutations type = gene set (sample n<7) or phenotype (sample n>7). Putative protein-protein interactions were identified using STRING [45]. Interferon-inducible genes were identified using Interferome v2.0 [46] Non-weighted Venn diagrams were generated using Venny 2.1 [47]. Hamster transcriptomic datasets were from 41. Comparisons were made based on ENSEMBL ids and gene symbol and annotation was used to compare hamster / mouse orthologs. Human transcriptomic data sets were from 39. Weighted Venn diagrams were generated in R using the eulerr package ( https://cran.r-project.org/package=eulerr) and correlation plots were created using the ggplot2 package ( https://cran.r-project.org/web/packages/ggplot2/index.html).

Results

Splenic response to L. donovani infection in BALB/c mice

As part of a larger program of research aimed at evaluating multiple aspects of EVL (the CRACK IT VIDR project; see Discussion), we infected two cohorts of female BALB/c mice with 2×10 7 amastigotes of L. donovani (CRACK-IT_1 and CRACK-IT_2) and at defined time points corresponding to increased parasite load and architectural changes in the splenic microenvironment ( Figure 1), spleen tissues were processed for genome-wide transcriptomic profiling by microarray. An overview of the expression data is presented as volcano plots and heat maps for naïve mice (total n=15) and infected mice (total n=25) examined at each time post infection (p.i.) ( Figure 2A and B). Some minor changes were noted in groups of naïve mice examined at each time point, despite confirmation of comparability in age, sex, weight and health status. To account for these minor variations, DE genes were identified by reference to time- and experiment-matched naïve controls ( Figure 2B). Using a mean FDR adjusted p value cut-off of 0.05 and a FC of 2 (Log 2FC = 1) for at least one of the infection time points, 9634 probes were scored as DE across the entire time series ( Table S1). The number of DE probes increased over time, peaking at d36 p.i. (d15: 1997 probes, Log 2FC -3.52 to +5.82; d17: 2826 probes, log 2FC -4.01 to +6.37; d21: 5124 probes, Log 2FC -4.54 to +7.59; d36: 5051 probes, Log 2FC -6.45 to +8.36; d42: 4946 probes, Log 2FC -6.30 to +7.13).

Figure 1.

Splenic response to L. donovani infection.

Spleens from BALB/c mice infected with L. donovani were removed at 15-, 17-, 21-, 36- and 42-days post infection, along with spleens from naïve age matched controls. A. Parasite burdens were determined by impression smears and are represented as mean LDU ± SE (n=5 mice per time point). B. Spleen weight for naïve and infected mice (mg ± SE). C– H. Representative histology of spleens from naïve mice ( C and D) and from d17 ( E and F) and d42 ( G and H) infected mice. Hematoxylin and Eosin x20 ( C, E and G) and x40 ( D, F and H); scale bars 50 µm. Blue arrows highlight parasite clusters. White dashed lines indicate red pulp (rp) / white pulp (wp) boundary. G and H show increased frequency of megakaryocytes, indicative of extramedullary haematopoiesis and loss of white pulp / red pulp discrimination. Data are pooled from two independent experiments covering the early and late phase of infection.

Figure 2.

Transcriptomic profiling of spleen response to L. donovani infection in BALB/c mice.

A. Volcano plots showing DE probes at each time point of infection (relative to matched naïve control mice). Data are shown as Log2 FC in expression against Log10 adjusted p value. B. Heat map of probe expression intensity across time series. Data are pooled for n=5 mice per group for clarity. Vertical blue trace indicates mean intensity signal. Sample “naïve 0” was used to calculate DE genes for days 15, 17 and 21 p.i. Samples “naïve 36” and “naïve 42” were used to calculate DE genes for their respective time points p.i.

After removal of multiple probes, 5096 annotated genes were identified within this DE probe set across all time points (d15, 1078 genes; d17, 1653 genes; d21, 2880 genes; d36, 2708; d42, 2805). To determine if genes related to specific biological responses were differentially regulated over time of infection, we used GSEA [34]. As shown in Figure 3A, Table 1 and Table S2, the most commonly enriched gene sets across the time course were interferon gamma signalling, TNF signalling, myogenesis, allograft rejection, IL-6 – JAK/STAT signalling, angiogenesis, G2M checkpoint, inflammatory response, oestrogen response, epithelial-mesenchymal transition, E2F targets, KRAS signalling and complement. Notably, the breadth of enriched pathways evolved over time. For example, pathways representing various metabolic processes, coagulation, apoptosis and hypoxia are more highly enriched late in infection ( Figure 3A). These changes occurred concurrently with loss of red pulp/white pulp differentiation and the onset of extramedullary haematopoiesis (as evident from the abundance of splenic megakaryocytes; Figure 1G and H).

Figure 3.

Characteristics of the spleen transcriptomic response to L. donovani infection in BALB/c mice.

A. Data were subjected to GSEA with reference to MSigDB Hallmark gene sets and normalized enrichment score over time is shown for the top 25 gene sets at any given time point. Days post infection are indicated by coloured dots (d15, black; d17, blue; d21, magenta; d36, brown, d42 red). B. IPA-predicted upstream regulators presented as a heat map based on z scores ranging from -9 (negative regulators) to +9 (positive regulators). C. IPA mechanistic pathway analysis indicating the 16 linked regulatory proteins, that together are predicted to regulate ~14% of all DE genes in the L. donovani infected spleen. Predicted relationships are indicated by lines (orange, leading to activation; blue, leading to inhibition; yellow, inconsistent; grey, not predicted). Box shading represents degree of predicted activation (orange) or inhibition (blue). D. Venn diagram showing distribution of 1645 Type I and II Interferon-induced genes, as identified using Interferome v2.01. E. Predicted Trim24 targets identified in IPA that show upregulated expression during infection. Predicted relationships are indicated by lines: Orange represents leading to activation; grey represents not predicted. Gene symbol box shading reflects extent of upregulation (red) and downregulation (green) in dataset. For an explanation of molecule shapes and relationship types, see http://qiagen.force.com/KnowledgeBase/articles/Basic_Technical_Q_A/Legend.

Table 1.

Summary of enriched gene sets in the spleen of L. donovani-infected BALB/c mice.

Data was subjected to GSEA against the MSigDB Hallmark gene set and represents a summary of the top 10 gene sets (by FDR) for each time point. For full dataset including gene lists, see S2Table. Normalised enrichment scores for all gene sets listed for all time points are shown in Figure 3A.	day 15	day 17	day 21	day 36	day 42
Gene Set Name	FDR q-value
HALLMARK_KRAS_SIGNALING_UP	7.80E-09	3.07E-11		2.70E-24	1.63E-20
HALLMARK_KRAS_SIGNALING_DN	2.67E-08
HALLMARK_ALLOGRAFT_REJECTION	1.15E-07	8.95E-10		9.97E-27
HALLMARK_IL6_JAK_STAT3_SIGNALING	2.48E-07
HALLMARK_INTERFERON_GAMMA_RESPONSE	2.40E-06	7.24E-26	1.07E-22	2.75E-40	1.88E-27
HALLMARK_ESTROGEN_RESPONSE_EARLY	9.16E-06
HALLMARK_INFLAMMATORY_RESPONSE	9.16E-06	3.07E-11	1.13E-19	4.08E-25	8.69E-23
HALLMARK_ANGIOGENESIS	1.06E-05
HALLMARK_COAGULATION	1.77E-05			1.90E-23
HALLMARK_EPITHELIAL_MESENCHYMAL_TRANSITION	3.20E-05	4.76E-13	1.13E-19	1.52E-18	3.09E-21
HALLMARK_INTERFERON_ALPHA_RESPONSE		1.26E-12
HALLMARK_HYPOXIA		1.26E-12		3.81E-20	3.38E-25
HALLMARK_MYOGENESIS		1.26E-12
HALLMARK_G2M_CHECKPOINT		3.07E-11	1.50E-36		1.96E-24
HALLMARK_TNFA_SIGNALING_VIA_NFKB		4.32E-09		2.37E-30	6.38E-26
HALLMARK_E2F_TARGETS			1.50E-36		1.63E-20
HALLMARK_MITOTIC_SPINDLE			1.13E-19
HALLMARK_IL2_STAT5_SIGNALING			6.16E-19	1.18E-22
HALLMARK_UV_RESPONSE_DN			2.57E-17
HALLMARK_XENOBIOTIC_METABOLISM			1.54E-14
HALLMARK_COMPLEMENT				4.08E-25	9.69E-20
HALLMARK_ESTROGEN_RESPONSE_LATE			1.54E-14		5.18E-22

We next used pathway analysis to look in more detail at the evolution of and the main upstream regulators of the DE genes identified at each time point post infection. Comparative analysis within Ingenuity Pathway Analysis (IPA) indicated that IFN γ, LPS, IL-1β, TLR4, TNF and STAT1 were amongst the highest scoring and significant upstream regulators predicted to be activated during infection (with activation z scores ranging from 2.89 to 9.5 and p values from ~10 -5 to 10 -48), whereas IL-10RA and TRIM24 were upstream regulators that were predicted to be down regulated (z scores of ~-4 to -9 and p values from 10 -7 to 10 -35; Figure 3B). It should be noted that predictions made by IPA are based on the expression patterns of genes that are known to be regulated by different upstream regulators and do not require that the upstream regulator itself appears as DE in the data set. For example, at day 36 p.i., activation of IFNγ was predicted by the changes in expression seen in 301 known IFNγ target genes, consistent in this case with an observed Log2FC change in expression of IFNγ of 3.474. IFNγ was further predicted to interact with 694 genes (or 13.7%) of all DE genes via a mechanistic network involving 15 other regulators ( Figure 3C). Similarly, using the Interferome database which contains expression data related to ~4000 interferon regulated genes (v2.01; [46]), 32.3% (1645/5096) of all DE genes scored as interferon inducible, of which 743 and 431 scored as uniquely regulated by Type I and II (IFNγ) interferons, respectively ( Figure 3D). In contrast, IL-10RA and TRIM24 were predicted to be inhibited based on the expression change of 108 and 41 genes, respectively, although in this case, neither of these predicted upstream regulators was itself DE. The 41 genes that predict inactivation of TRIM24 during L. donovani infection are shown graphically according to subcellular localisation in Figure 3E,. The lack of differential expression of IL-10RA and TRIM24 may reflect that function is not transcriptional regulated or that transcriptional regulation was below the Log2FC cutoff employed in this analysis (due to low abundance of cells expressing the gene of interest or only minor changes to mRNA abundance).

Finally, we examined expression of cytokine and chemokine genes as well as genes for their respective receptors, to provide a high-level view of immune pathways activated during infection in this target tissue ( Figure 4). Cytokines with increased mRNA accumulation at all time points included Il10 and Il21, whereas for others (e.g. Il1b and Il6) there was progressive increase in mRNA accumulation over time. Ifnγ mRNA abundance also increased over time post infection, whilst enhanced accumulation of Tnf mRNA was observed only late in infection and other TNF superfamily members showed variable responses. Distinct temporal regulation of chemokine mRNA abundance was also observed. In contrast to chemokine ligands, mRNAs for many chemokine receptors were found at lower abundance in infected mice compared to controls (e.g. Ccr1 and Ccr9), and the same trend was also observed for some cytokine receptors, notably Il17rb, Il17rd and Il17re.

Figure 4.

Cytokine and chemokine ligand and receptor gene expression in the spleen of mice infected with L. donovani.

A– C. Heat maps representing Log2 fold change in mRNA abundance for chemokines, interleukins and IFNγ / TNF family members and their receptors are shown. Only genes that were DE for at least one time point post infection are shown. Blue blocks represent down-regulated genes and red blocks indicate up-regulated genes. White blocks represent genes that were not significantly DE at a given time point.

Comparative analysis of splenic response in mice and hamsters infected with L. donovani

In contrast to murine EVL, hamsters infected with L. donovani develop progressive and fatal visceral leishmaniasis, providing a better model of end stage human disease. The heightened sensitivity of hamsters to L. donovani infection has been attributed, in part, to a deficiency in NOS2 production, resulting from a promoter polymorphism similar to that seen in humans [48]. In addition, alterations in immune response associated with T cell exhaustion and myeloid cell dysregulation have been reported based on transcriptomic as well as functional analysis [41, 42]. However, a detailed comparison of the pathophysiology in infected hamsters and mice has not been conducted. We made use of recently published RNA-Seq data derived from the spleens of hamsters infected with L. donovani infection [41] to compare these two infection models, albeit with the caveat that different technologies were used for measuring mRNA abundance. Kong et al. [41] observed 4360 DE transcripts at d28 post infection, of which 2340 (53.7%) were up-regulated and 2020 (46.3%) were down-regulated compared to spleens from naïve hamsters. In comparison, in L. donovani infected mice at the peak of splenic infection in our study (d36 p.i.), we identified 41.2% of DE probes as upregulated and 58.8% as down-regulated ( Table S1). By filtering on FDR and FC (<0.05; >2FC), removing the bottom quartile of hamster genes that were expressed with a cpm of <1 and by averaging across multiple probes/reads, we generated three lists of genes, those DE in hamster and not mouse (1504 genes), DE in mouse and not hamster (2239 genes) and those DE in both (485 genes, Figure 5A and Table S3).

Figure 5.

Comparison between DE genes in hamster and mouse spleen.

DE genes identified in mouse spleen (d36 p.i.; this study) were compared to those identified in 41. A. Venn diagram showing overlap of DE genes. B. Correlation plot of Log2FC for mouse and hamster DE genes. C– F. STRING analysis of predicted protein interactions for DE genes identified as upregulated in both mouse and hamster ( C; yellow = cell cycle, blue = immune response), down-regulated in mouse and hamster ( D; red = anatomical structure, blue = vascular development), up-regulated in hamster only ( E; yellow = negative regulation, blue = immune system), and down-regulated in hamster only ( F; green = regulation of multicellular processes, red = complex of collagen fibres, yellow = drug metabolism, blue = serine proteases trypsin domain). All DE genes are shown on the images to provide a visual representation of the extent to which these genes are clustered into common pathways. For full details of all DE genes, see Table S3.

Amongst common DE genes, there were those that were up or down consistently in both species, and in more limited number, some that showed differential expression in alternate directions (35 up in mouse and down in hamster and 37 up in hamster and down in mouse; Figure 5B and Table S3). To further examine these DE gene lists, we used STRING, a protein-protein interaction database [45]. We identified no discernible functional pathways associated with the small number of inversely DE genes. Analysis of common upregulated DE genes identified a highly significant enrichment for immune system process (GO:002376; FDR 7.37×10 -23) and cell cycle process (GO:0022402; FDR 2.75×10 -19) related genes that form two major clusters amongst all commonly upregulated genes ( Figure 5C). Amongst common down-regulated DE genes, although there was no major clustering of genes there was a significant enrichment in GO terms related to changes in morphology, pathways for anatomical structure morphogenesis (GO:0009653; FDR 3.22×10 -14) and vascular development (GO:0001944; 1.82×10 -11), the latter including negative regulators of angiogenesis such as Mmrn2 and promoters of angiogenesis such as Vegfa ( Figure 5D). These data are consistent with the notion that the well documented changes in vasculature associated with murine infection [9] also occur during hamster infection.

Amongst DE genes upregulated only in infected hamsters, GO pathways related to negative regulation (e.g. GO: 0048519; FDR 1.13×10 -7) and immune system processes (e.g. GO:0002376; FDR 2.1×10 -7) predominated, again forming significant clustering in the dataset ( Figure 5E). Most notable amongst the top 25 up regulated genes were Ido1, a potent negative regulator of immune responses (245-fold upregulated) and Arg1, associated with Th2-associated cytokines (36-fold upregulated). Surprisingly, Il10 was not upregulated in splenic tissue. In contrast, DE genes downregulated in hamsters only related to morphogenesis and organismal processes (e.g. GO:0051239: FDR 5.78×10 -9), complexes of collagen trimers (GO:0098644; FDR 3.4×10 -11), drug metabolism (KEGG:00983; FDR 7.94×10 -8) and serine proteases, trypsin domain (IPR001254; FDR 0.00741). Of the top 25 down-regulated genes, 10 encoded amylase or amylase precursor proteins and 10 have serine peptidase activity ( Table S3).

Amongst DE genes upregulated only in murine infection, genes associated with GO terms related to host defence and immune response again predominated (e.g. GO:0006955; FDR 2.03×10 -10; GO:0006952; FDR 4.96×10 -10). Notable amongst the top 25 upregulated genes were Ptgs2 (or Cox-2), a key enzyme in the pathway leading to prostaglandin synthesis (124-fold upregulated), Tff1, associated with chronic inflammation (75-fold upregulated), Ctsg, thought to promote excessive inflammation (33-fold upregulated) and the neutrophil related proteins Lcn2 and Ngp (19.7 and 24.3-fold upregulated respectively). Notable in the top 25 downregulated genes were Cd27 (langerin; 23-fold downregulated), the sodium voltage gated channel encoding Scn3a (18-fold down regulated), Cd209c (SIGNR2; 9-fold), and Il22ra2, a soluble receptor that inhibits IL-22 activity (12.7-fold).

Finally, given its importance to the host response to infection, we examined the status of macrophage activation, and observed a good degree of concordance for M1-associated transcripts (as described in Figure 6 of 41) between hamster and mouse (55% of genes, including Ccl2, Cxcl9, Cxcl10, Ifng, Irf1, Irf7, Irg1, Socs3, Stat1 and Ccr7). In contrast, M2-associated genes were less commonly regulated in mouse compared to hamster spleen (23% of genes listed in Figure 7 of 41, including Ccl2, Cd209a, Cebpb, Chi3l1, and Il1rn). These data are in accord with the above analysis indicating an enhanced Th-2 type inflammation in hamster compared to murine EVL.

Hepatic response to L. donovani infection in BALB/c mice

We adopted a similar approach to examine transcriptional changes in the liver of L. donovani infected mice, where the process of granulomatous inflammation has been well-characterised at the histological level [13, 49, 50]. In contrast to the persistent parasite burden observed in the chronically infected spleen, parasite burden was significantly reduced in the liver by d42 p.i. ( Figure 6A). This was commensurate with the development and maturation of host protective granulomas, as indicated by quantitative morphometric analysis of granuloma number and area ( Figure 6B–F) and through observation of H&E stained sections ( Figure 6G–J).

Figure 6.

Hepatic response to L. donovani infection.

Livers from BALB/c mice infected with L. donovani were removed at 15-, 17-, 21-, 36- and 42-days post infection. A. Parasite burdens were determined from impression smears and are represented as mean LDUs ± SEM. B. Granuloma density (mm 2; green dots) and average granuloma size (µm 2; yellow squares) are shown per mouse and with group mean. C– F. Distribution plot of granuloma size (ranging from 150 µm 2 to 6140 µm 2) for 15, 21, 36 and 42 days p.i. respectively. Inset shows representative granulomas identified using TissueQuest, with F4/80 staining (red) and attributed size in µm 2. G– J. Granuloma morphology at 15, 21, 36 and 42 days p.i. respectively. H&E x40 magnification; scale bars 50 µm. Blue arrows indicate amastigote clusters. Dashed line denotes granuloma perimeter.

In the infected liver, 16192 probes were scored as DE (>2 FC, 5% FDR) for at least one time point ( Figure 7 and Table S4). As seen in volcano plots ( Figure 7A vs Figure 2A) there was a tendency towards greater fold changes in mRNA accumulation compared to spleen with a bias towards up-regulated genes.

Figure 7.

Transcriptomic profiling of hepatic response to L. donovani infection in BALB/c mice.

A. Volcano plots showing DE probes at each time point of infection (relative to matched naïve control mice). Data are shown as Log2 FC in expression against Log10 adjusted p value. B. Heat map of probe expression intensity across time series. Data are pooled for n=5 mice per group for clarity. Vertical blue trace indicates mean intensity signal. Sample “naïve 0” was used to calculate DE genes for days 15, 17 and 21 p.i. Samples “naïve 36” and “naïve 42” were used to calculate DE genes for their respective time points p.i.

By GSEA, we identified 20–24 hallmark gene sets that were associated with infection over the time course of infection (using a conservative threshold of 5% FDR); these enriched gene sets showed a high level of concordance and included allograft rejection, inflammatory response, IFNγ and IFNα signalling, TNF- , IL-2- and IL-6-related signalling pathways, pathways related to cell cycle, apoptosis, and complement Several were also found associated with the response to infection in the spleen ( Table 1 and Figure 8). There was, however, liver-specific enrichment for hallmark gene sets representing regulation of Kras signalling, epithelial and mesenchyme transition, angiogenesis, notch signalling and apical junction and apical surface. The latter are likely to reflect underlying changes in hepatocytes, which although not grossly abnormal in structure in H&E sections, may nevertheless have an important role in liver pathophysiology during disease [51].

Figure 8.

Comparison between hepatic and splenic responses to L. donovani by GSEA.

Data were subjected to GSEA with reference to MSigDB Hallmark gene sets and normalized enrichment score (NES) over time is shown as a heat map for significant gene sets enriched in liver spleen and blood at the state time points.

mRNA accumulation for genes associated with cytokines, chemokines and their associated receptors were also compared over time ( Figure 9). In general, the breadth and intensity of response seen in the liver was greater than in the spleen, a result at least in part due to the greater change in cellular composition between resting and inflamed liver compared to resting and inflamed spleen. Nevertheless, it is also likely that some of the cytokine / cytokine receptor pathways that are more pronounced in the liver may reflect the self-cure phenotype associated with this organ. For example, at least at the level of mRNA abundance, IFNγ and IL-10 are not discriminatory between spleen and liver, whereas early increases in mRNA abundance for TNF and a broad range of TNF superfamily members appears liver specific. Of note, for many cytokines and chemokines, day 36 represented the peak of the response in the liver, with reduced mRNA accumulation at day 42 associated with the reduction in parasite load seen at this time ( Figure 6).

Figure 9.

Cytokine and chemokine ligand and receptor gene expression in the liver of mice infected with L. donovani.

A– C. Heat maps representing Log2 fold change in mRNA abundance for chemokines, interleukins and IFNγ / TNF family members and their receptors are shown. Only genes that were DE for at least one time point post infection are shown. Blue blocks represent down-regulated genes and red blocks indicate up-regulated genes. White blocks represent genes that were not significantly DE at a given time point.

Finally, we examined the expression of inhibitory receptors associated with T cell exhaustion [5, 52, 53], to determine whether expression at the whole tissue level provided any clues as to the rate of cure for L. donovani infection seen in these two organs ( Figure 10). Inhibitory receptors were more prominently associated with the hepatic response, again likely reflecting the more dramatic increase in relative T cell frequency in liver vs. spleen. In spleen, Lag3 was the only inhibitory receptor that remained elevated over the entire course of infection. Of note, similar to mouse spleen, Lag3 was also significantly upregulated in hamster spleen whereas Havcr2 (Tim-3) was also significantly down-regulated ( Table S3).

Figure 10.

Expression of immune checkpoint genes in the spleen and liver of mice infected with L. donovani.

Heat map representing Log2 fold change in mRNA abundance for immune checkpoint genes in spleen and liver of L. donovani-infected BALB/c mice. Only genes that were DE for at least one time point and at least one organ are shown. Blue blocks represent down-regulated genes and red blocks indicate up-regulated genes. White blocks represent genes that were not significantly DE at a given time point.

Whole blood transcriptional response to L. donovani infection.

As whole blood is the most accessible target tissue for transcriptomic profiling in humans and previous work had suggested that whole blood signatures might be useful avatars of the systemic response, we compared the response between spleen, liver and blood. 26,836 probes (12348 annotated genes) were identified as DE (FDR 0.05, FC2) in at least one time point ( Table S5). By GSEA blood was, not surprisingly, more similar to spleen than liver ( Figure 8). We next took advantage of a recently published whole blood transcriptomic analysis [39] to compare responses across murine EVL and human VL. Using day 36-infected mice as a comparator, our results suggest a high divergence in whole blood response between EVL and human VL ( Figure 11).

Figure 11.

Comparison of whole blood transcriptome in murine EVL and human VL.

A and B. Murine whole blood data from day 36 p.i. was re-analysed with a cut off of log2FC= 1.5 to make comparable with the human data presented in 39. Venn diagram indicates a minimal overlap between DE genes that are up regulated ( A) and down-regulated ( B) during infection. C and D. STRING analysis of the common up-regulated ( C) and down-regulated ( D) genes. E and F. Most significant (by FDR) GO terms and Kegg pathways associated with mouse and human specific up-regulated ( E) and down-regulated ( F) genes, generated by STRING.

Identification of common response signature in spleen, liver and blood during murine L. donovani infection

Finally, we asked whether it was possible to identify a common response signature found in all tissues that might act as a biomarker of infection and which, when measured in blood might reflect or predict systemic disease. When analysed for genes that were DE at all time points in all tissues, we did not identify any genes that were significantly down regulated in expression. In contrast, 26 genes were consistently and significantly up-regulated across blood, spleen and liver at all times p.i in murine EVL ( Figure 12). These genes, of which 25/26 (96%) are interferon responsive, clustered in STRING analysis around 2 nodes ( Cxcl9 and Gbp1: Figure 12B and C) and reflected a significant enrichment for GO TERMS including “cellular response to interferon-beta” (GO:0035458, FDR 7.19×10 -8), “defence response to protozoans” (GO:0042832, FDR 7.22×10 -8), “immune response” (GO:0006955, FDR 9.71×10 -8), and symbiont containing vacuole” (GO:0020003, FDR 2.15×10 -9). The gene set also contained 3 apolipoprotein L genes ( Apol10a, Apol10b and Apol11b) which are phylogenetically related amongst the 12 members of the murine Apolipoprotein L family [54]. Of note, 19/26 (73%) genes in this signature were also upregulated in hamster spleen ( Figure 12B), whereas 6/26 (23%) were upregulated in human whole blood [39]. Common to all data sets were genes encoding IFNγ, the CD20-like antigen membrane spanning 4-domains A6A (MS4A6A/ Ms4a6d), and the interferon-inducible genes guanylate binding protein 2 (GBP2/ Gbp2), suppressor of cytokine signalling (SOCS1/ Socs1) and tryptophanyl tRNA synthetase (WARS/ Wars).

Figure 12.

A 26-gene transcriptional signature is present in blood, spleen and liver throughout the course of L. donovani infection in BALB/c mice.

A. Venn diagram showing relationship between genes consistently up-regulated across time in blood, spleen and liver. B. Heat map showing expression over time for the 26 gene signature in each tissue. Data are shown as Log2FC, indicated by the scale bar. Asterisks denote that genes were also upregulated in the spleen of hamsters infected for 28 days (data from 41.) C. STRING analysis of the 26-gene signature indicates two hub genes ( Cxcl9 and Gbp1) forming a network that encompasses 21/26 (81%) of these genes.

Discussion

The immune response to L. donovani infection has been extensively analysed using a wide array of technologies, including use of gene KO mice, advanced imaging techniques and cell/cytokine manipulations [5, 8, 50, 55– 57]. Yet a comprehensive picture of the dynamic complexity of the host response has not previously been reported. Here, we provide data from a detailed transcriptional analysis of the response of BALB/c mice to L. donovani infection and compare this where possible to existing data sets from other species, including humans. This high-level overview of transcriptional changes associated with infection poses a number of important questions and provides an important additional data resource for the leishmaniasis research community.

BALB/c mice were chosen for this study due to their extensive use as an immunological model of infection and as a tool for the pre-clinical evaluation of drugs and vaccines [58– 61]. In this and most other mouse strains on a Slc11a1 (Nramp1) mutant background, such as C57BL/6, hepatic infection is self-resolving and associated with the development of granulomatous inflammation, whereas splenic parasite load persists, often for several months, accompanied by marked alterations to organ microarchitecture [13]. We selected time points when hepatic and splenic parasite burden were increasing and also to reflect the time frame over which a dichotomy in host resistance becomes apparent in spleen and liver. Where possible, we have also mined existing data to look for similarities and/or differences in transcriptional response. Whilst these comparisons are informative, some degree of caution needs to be exerted on their interpretation. Platforms used for analysis of transcriptomic changes differ (Agilent gene arrays, Illumina Bead Chip arrays, RNA-Seq), with microarrays often having more limited dynamic range than RNA-Seq. For this reason, we have not attempted any quantitative comparison in expression levels across platforms. Rather, our comparisons are drawn on the basis of comparing DE gene lists generated using standard analysis and statistical tools suited for each platform. In addition, variables such as tissue parasite load and time of infection are difficult to standardise (particularly for human disease where the time of infection is usually unknown), and with the exception of our study, data is currently restricted in most species to a single tissue site. It should also be remembered that whether comparing between tissues in the same host or across species, lists of DE genes may reflect changes due to either alterations in gene expression pattern of individual cells or due to changes in tissue cell composition. Although in some cases our results were unexpected and showed dissimilarity between experimental models and human disease, this should not be taken as an argument against the value of research in such models. Rather, it stresses that once proof-of-concept has been achieved in a model, additional validation based on human data is essential. Furthermore, such human data may form the basis for selecting an appropriate pre-clinical model for the question at hand, or indeed for engineering more refined models of novel aspects of human disease.

Splenic histopathology of infected BALB/c mice, as expected, indicated extensive remodelling of the splenic architecture and the late onset of extramedullary haematopoiesis, similar to that seen in other mouse strains [8]. Coincident with these changes, we observed transcriptional signatures associated with the activation of immune responses, angiogenesis and coagulation pathways. Interferon-related pathways were highly dominant. Approximately 13% of the identified DE genes were predicted to be downstream of IFNγ signalling and ~30% of all genes were scored as interferon inducible using the Interferome database [46]. IPA identified two candidate negative regulators IL-10Ra and TRIM24. IL-10 signalling through IL-10Ra (CD210) is well recognised for its role in inhibiting immune clearance of Leishmania in mouse models and in humans [1, 2, 7, 62– 64], operating to directly regulate macrophage activation and also the function and development of IL-10 producing Th1 and Tr1 cells [1, 63– 66]. Therapeutic interventions to target IL-10 mediated immune regulation in human VL have been posited [67]. TRIM24 is a bromodomain-containing transcriptional regulator protein that has recently been the focus of attention as a regulator of STAT3 recruitment during oncogenesis [68] and as a positive regulator of Th2 cytokine production in Th2 cells [69]. The predicted inhibition of TRIM24 would, therefore, be consistent with a dominant Th1 type immune response during murine EVL. TRIM24 is also observed as a target for inactivation during infection with lab-adapted H1N1influenza A virus [70] and has been shown to be targetable using novel heterobifunctional protein degraders [71], suggesting that pathways controlled by TRIM24 could be further manipulated in favour of host protection.

In contrast to the murine EVL model, data generated by Melby and colleagues suggests that the immune response in hamster EVL favours a more prominent Th2 like response [40, 41]. In addition, their analysis also indicates that macrophages in hamster EVL are conditioned to respond to activation with the induction of a program of gene expression that favours amastigote growth, including Arg1, Ido-1 and Irg1 and associated with STAT3 activation [41]. At the peak of splenic infection in our study (d36 p.i.), which we considered the best approximation of the hamster model that we could achieve based on the degree of pathology, we did not observe increased accumulation of Ido1 in murine EVL and a comparison of signature M1 and M2 genes suggested that in murine EVL, M1 activation was favoured. At a tissue level, alternate regulators of inflammation such as Ptsg2 and Il-10 appeared to be favoured in murine EVL. Upregulation of Ptsg2 (Cox2) is consistent with the enhanced production of prostaglandins in spleen cells from mice with L. donovani infection [72] and with the host protective effect of Cox2 inhibition in vivo [73].

There has been considerable recent interest in the potential for host-directed therapies in different forms of leishmaniasis [4, 5, 67, 74– 77]. Of interest, therefore, is how well pre-clinical models may predict novel approaches applicable to human VL. Our analysis showing increased mRNA abundance for a variety of checkpoint inhibitors including CTLA4, PD1, LAG3 and others ( Figure 10) supports the notion that these may be candidates for the development of immunotherapy, along with key suppressive cytokines like IL-10. In the mouse, these checkpoint inhibitors (along with many cytokines, chemokines and their receptors) have a complex pattern of temporal and tissue-specific regulation which, together with the potential for cell migration between sites, makes predictions concerning the impact of immunotherapy challenging. Some of these checkpoint targets are equally represented in mouse and hamster spleen and human whole blood e.g. LAG3 and CTLA4, whereas others show species restricted regulation. For example, IDO-1 is up regulated in hamster spleen and human whole blood, but not in any mouse tissue examined. Importantly, these diverse data as well as the clear dichotomy of expression at a tissue specific level, indicate the need to build a more comprehensive and uniform knowledge base regarding these pathways in different models and forms of human disease.

Analysis of whole blood has been a valuable tool for monitoring transcriptional changes associated with infection and even for predicting drug efficacy [19], but the relationship between transcriptomic changes in blood and tissue remains poorly understood. Pooling data across all tissues and time points, we were able to identify a 26-gene signature that was common to spleen, liver and blood during murine VL. 96% of these genes have been reported as interferon inducible, clearly indicating the dominant nature of interferon signalling during murine VL and that this common signature reflects disease in all organs/tissues studied. The value of defining signatures of this type might be to understand common pathways underlying inflammation or to provide new biomarkers that reflect the level of ongoing systemic disease that can be monitored in blood. We will further elaborate on the use of this signature to evaluate chemotherapy induced cure in murine VL in a subsequent manuscript (Forrester et al. manuscript in preparation). The clinical value of this signature is however questionable at this stage, given the lack of similarity between whole blood transcriptomic responses observed in our study and that of humans with VL in Brazil [39]. In a similar study in acute melioidosis, 26.9% of genes were similarly regulated in murine and human infection [21]. The reasons for this lack of similarity in VL are likely to be multiple, including the chronicity of infection, difference in pathogen ( L. donovani vs. L. infantum), specific characteristics of Brazilian VL compared to VL on other continents [12] or differing analytical methods. It would be informative, given the above discussion, to compare human whole blood from a greater range of patient populations with VL to better understand the extent of temporal and geographic diversity in the transcriptome of patients with VL. Likewise, data on hamster whole blood might help ascertain whether the blood transcriptome in this model more closely resembles that found during human disease.

Finally, a deeper understanding of subtleties of the immunopathology of leishmaniasis in different target tissues (e.g. granulomatous vs. non-granulomatous) and their relationship to gene transcriptional changes will require more in-depth analysis than that reported here, using new techniques to integrate image analysis with complex “omics” and other meta data. Such approaches are fuelling significant advances in precision medicine in other fields, notably cancer and neuroscience [78– 80], but have to date had limited application in the field of infectious disease. To facilitate the development of such approaches in leishmaniasis, we have generated a digital whole slide collection of the tissue sections generated in this study, stained with both H&E and markers to identify myeloid cells (F4/80 and 1A8). These are available via a new network designed to support the development of digital pathology in leishmaniasis ( www.leishpathnet.org), and improve transparency and data sharing. In addition, the data collected from the animals used in this study forms part of a larger program of work (the CRACK IT Virtual Infectious Diseases Research Challenge; https://crackit.org.uk/challenge-16-virtual-infectious-disease-research) that had the goal of generating an in silico model of the immune response during murine EVL, with the aim of reducing the number of animals required for basic research on immune regulation and for evaluating potentially novel therapeutic interventions ( www.leishsim.org). Future manuscripts from the CRACK IT program will extend the transcriptomic analysis presented here to include comparative studies with L. infantum infection, the response of L. donovani-infected BALB/c mice to chemotherapy, and the response to L. donovani infection in C57BL/6 (B6) and B6. Rag1 -/- mice.

Data availability

Microarray data is available via the Gene Expression Omnibus (Ascension number GSE113376: https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE113376).

OSF: CRACKIT Virtual Infectious Diseases project. Supplementary Table 1– Supplementary Table 5 (see information below in Supplementary material) and raw data for Figure 1A and B and Figure 6A–F are available, https://doi.org/10.17605/OSF.IO/9WSDK [81].

Data on OSF are available under the terms of the Creative Commons Zero "No rights reserved" data waiver (CC0 1.0 Public domain dedication).

Whole slide images and individual mouse metadata are available from www.leishpathnet.org (study designations CRACKIT-1 and CRACKIT-2). Requests for access to tissue samples from these studies will be accommodated where possible and subject to availability.

RNA-Seq data from a hamster model of VL is available from: https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE91187

Blood transcription of human infection are available from: https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE77528

The authors satisfactorily dealt with all the questions placed in our previous analysis.

We have no further comments to make.

We have read this submission. We believe that we have an appropriate level of expertise to confirm that it is of an acceptable scientific standard.

This study provides an analysis of transcriptional changes occurring in the Leishmania donovani infected mouse over time in the spleen and liver. There is no comparable kinetic analysis of transcriptional changes in visceral leishmaniasis, thus indicating the significance of the work. The authors also attempt to compare their findings with published transcriptional analyses in hamster and human VL, which is useful. Overall, the study is well-done with some important caveats that the authors are aware of and which they mention in the text. First and foremost, the kinetic study in mice is a combination of two experiments and the authors might be clearer in how they handled the potential batch effects that one might anticipate seeing. Looking at Figure 2 it would appear that there are significant differences between the first experiment (Day 15, 17 and 21) and the second experiment (day 36 and 42).  Second, a comparison between data generated by microarray and RNAseq will have several notable issues, and while the authors mention this issue they might expand a bit on how those differences will influence interpretation of their results. It can be useful to make these comparisons, but critical to ensure the readers understand the significant caveats. Finally, it would be helpful to have a clearer discussion of what they have found out about tissue differences in response to infection between the spleen and liver, and how those differences might be responsible for differences in the parasite control that is seen in these two tissues.

The minor problem that becomes a bigger problem as it is evident in many of the figures is the inability to read several of the figures because of font size, as well as some figures that may not be informative for the readers. For example: Fig. 1C-D: arrows almost invisible: Fig. 3A. difficult to evaluate; Fig. 3C- Font size; Fig. 3E- Font size; Fig 5B (a really cool figure)- but font size too small for the labeled genes; Fig. 5C-F- difficult to interpret; Fig. 11C, D- Font size; Fig. 12C- Font size. One additional question: There appear to be no changes in the granulomas in Fig. 6B from day 21 to 42. Is this expected?

In summary, in spite of the caveats mentioned above, this study provides data that can be mined by others to better understand the immune responses in VL.

I have read this submission. I believe that I have an appropriate level of expertise to confirm that it is of an acceptable scientific standard.

1. batch effects resulting from conducting two independent experiments

We are aware of the possibility of batch effects and have tried to control for these as best as possible. The size and scale of the experiments performed for this study (and also to support additional analysis to be described elsewhere) made it logistically impossible to use a single cohort processed on the same day.  We confirmed health status across cohorts according to FELASA guidelines and had matched “naïve” control mice for each cohort (and for each time point in CRACKIT_2). Differential gene expression was conducted relative to these matched controls rather than to a common “naïve” group to further reduce batch effects. These details are described in the M&M and Results.

2. issues resulting from comparison of RNA-Seq to microarray

We have expanded our discussion on this point as suggested by the referee

3 tissue specific differences in parasite control in spleen and liver

We have included some additional comments on this issue throughout the manuscript.  However, as this data represented the platform on which to base further mechanistic studies, we have refrained from making too many conclusions purely on expression data.

4. font size of figures

We have reviewed all the figures for legibility and where necessary either increased the font size, of re-designed the figure to improve readability. For the STRING diagrams in Figs 11 and 12, the main aim is to show the extent to which different groups of DE genes are clustered. We tried increasing font size, but this actually makes the figure less clear due to proximity of each name. We think it valuable for readers to be able to view the gene names (even if requiring zoom on a pdf version of the manuscript and so have retained a small font on Figs 11 and 12).

5. was no difference in granulomas between d21 and d42 expected?

There are some subtle differences in the size distribution over time, and at d42 there is already some killing of parasites within granulomas (as evident from Fig 1A). This is not fully reflected in change in granuloma size, however.

The manuscript entitled “Tissue and host species-specific transcriptional changes in models of experimental visceral leishmaniasis” by Ashwin et al. describes tissue specific alteration in transcriptome in Leishmania donovani infected Balb/c mice in a longitudinal infection and compares the results with those from other species.

The number of DE probes was studied between normal and infected splenic tissues of Balb/c mice as a function of time, peaking at d36.  After removal of multiple probes, the number of annotated genes was 5096 set across all time points. The most commonly enriched genes were IFN-Y, TNF-α signalling, IL-6-JAK/STAT signaling, E2F targets, KRAS signaling and others. Interestingly these changes are coupled with extra medullary myelopoiesis in spleen. IPA showed significant upstream regulators predicted to be up regulated during infection whereas IL-10RA and TRIM24 were upstream regulators that were predicted to be down regulated, despite not being DE themselves. Using the Interferome database they could show that 32.3% of all DE genes scored as interferon inducible. Further pathway analysis was done to find upstream regulators of the DE genes. The mouse data was compared with the hamster and human data available in the literature. It is also highly appreciable that the data collected from the animals used in this study forms part of a CRACKIT Virtual Infectious Diseases Research Challenge program that focuses on generating an in silico model of the immune response during murine EVL in order to reduce the number of animals required for basic research on immune regulation and for evaluating potentially novel therapeutic interventions. Besides this, they have also made Microarray data accessible in the public domain through Gene Expression Omnibus a part of NCBI Databases.

And finally, they have confirmed that there is marked tissue-specific alteration in the transcriptome of infected mice over time and also identified commonality and differences between murine, hamster and human responses to infection that were previously unrecognized. They have shown the commonality of interferon-regulated genes while confirming a greater activation of type 2 immune pathways in infected hamsters compared to mice, and also pointed out that Cytokine genes and genes encoding immune checkpoints were markedly tissue-specific and dynamic in their expression and pathways. Through data, they have made an effort to address the value of measuring peripheral blood transcriptomics as a potential window into the underlying systemic disease. They could show a 26-gene signature that is common to spleen, liver and blood during murine VL. Oddly enough, there is a lack of similarity in whole blood transcriptomic response between mouse and human systems.

The major challenge now is the absence of an experimental model that would mimic the spectrum of human disease. Studies on the animal model provide us with a wealth of information about the host-parasite interaction, but there is no such model available that would reproduce human disease. Thus, so many questions remain unanswered in order to design control strategies. This observation re-enforces the notion that there is a gap between experimental data and human disease on leishmaniasis.

General comments:

The authors try to compare DE transcripts at d28 post infection from naïve hamsters to L. donovani to 36 days post infected mice and their study is based on the best match (page 9 of 21). It will be nice if the basis for the match criterion can be elaborated.

Figure 5A is a Venn diagram showing the overlap of DE genes in page 10. It looks like a weighted Venn diagram and the authors have cited Venny 2.1 for making Venn diagram but I wonder if the weighted Venn cannot be formed by Venny 2.1.  It will be helpful if the tool or package that has been utilized for its formation can be included. In case it is just a normal Venn diagram then the figure requires an edit.

Figure 5B: Please provide the citation to the tool or package that has been utilized for the correlation plot of Log2FC for mouse and hamster DE genes.

There are a few issues like parasite stain, duration of infection and resulting host-pathogen interaction etc that are beyond the authors’ control. This raises the question whether even monkeys can be used for testing vaccine candidates as little is known if the immune response to leishmania is similar to that observed in human.

Nevertheless as a whole I think it is a great effort by the group. It is a pleasure to read the article.

I have read this submission. I believe that I have an appropriate level of expertise to confirm that it is of an acceptable scientific standard.

1. elaborate on why mouse d36 and hamster d28 was considered “best match” for analysis

We have amended the text in the Results to indicate “at the peak of splenic infection” rather than “best matched” and made a further comment to this effect in the Discussion.

2. software to generate weighted Venn diagram in Figure 5A

We have added this information to the M&M

3. software to generate correlation plot in Figure 5B

We have added this information to the M&M

This is a very informative paper that compares gene expression during the early course of EVL in mice and hamsters, the differential expression pattern in spleen, liver and whole blood and the experimental data with the available data obtained from HVL. The strategy used is robust and very transparent. The results are in accordance with the methodology used. It will help to define new strategies for using of experimental models and to define new markers of immune response for HVL.

METHODS:

Why did the authors perform the study using two different experiments CRACK-IT_1 and CRACK-IT_2?

Tissue transcriptomics: what impact may the early amplification and removal of outliers by normalisation have had on the comparative study presented in the manuscript?

Was any additional normalization needed for comparisons performed between the data presented in this manuscript and those presented in references 39 and 41?

RESULTS:

Figures 1A and 1B_Raw (Supplementary data) The mean spleen weight of control animals at 36 days was ~4x higher than of animals with 14 and 42 days? Was the histology of these enlarged spleens similar to that observed in the other control animals? Do the authors think that including animals with enlarged spleen may have interfered with the data on gene expression?

Figure 1C-1H:

Technical details: For some reason the figures present some degree of opacity impairing adequate visualization. The blue arrows are so small that is difficult to find and the details such as infected macrophages are barely seen.

For the usual reader the loss of red pulp/white pulp differentiation mentioned by the authors in the text may not be easily seen. A more representative image and addition of markers and legends indicating these morphological changes may help to present the data.

Still about the figure comments presented in the text. Extramedullary splenic hematopoiesis is common in normal mice and has a tendency of increasing with age. Did the authors make any comparative quantitative assessment of spleen hematopoiesis between control and infected animals?

Figure 6B: please replace “green square” with “green dots”.

Figure 6E and 6F: I had the impression that there was a substantial decrease in F4/80+ cells in the granuloma. Was there a change in cell population in the granuloma at these stages of infection?

Figure 10: It is interesting that Ptger4, Klrag9 and Havcr2 are upregulated in liver and down regulated in spleen. Can the author comment on this data?

We have read this submission. We believe that we have an appropriate level of expertise to confirm that it is of an acceptable scientific standard, however we have significant reservations, as outlined above.

Methods

1.  why two different experiments for the time course analysis?

The size and scale of the experiments performed for this study (and also to support additional analysis to be described elsewhere) made it logistically impossible to use a single cohort infected and processed on the same day.  To minimise batch effects, we confirmed health standard across both cohorts according to FELASA guidelines and had matched “naïve” control mice for each cohort (and for each time point in CRACKIT_2).  Differential gene expression was conducted relative to these matched controls rather than to a common “naïve” group to further reduce batch effects. These details are described in the M&M and Results.

2. impact of early amplification and removal of outliers by normalisation on the comparative study

There are multiple reasons to be cautious in conducting comparative studies with data derived from different platforms and we discuss these in the text.  It is difficult to know the significance of these effects, but we have erred on the cautious side e.g. avoiding direct quantitative comparisons and commenting only on DE genes identified by conventional platform specific approaches.

3. was any additional normalisation used to compare mouse, hamster and human data sets?

We have added a further statement to the M&M to clarify the text already in the Results section.

Results

1. Fig1A/B Raw data on spleen size in naïve controls at d36

We thank the reviewer for pointing out this transcribing error in the raw data file. This has now been corrected (naïve spleen weight at d36 = 86+/-5 mg)

2. improve clarity of Figure 1C-H (arrows etc, demarcation of red and white pulp)

The figures have been corrected as suggested to improve clarity

3. comparison of spleen hematopoiesis between control and infected mice

We have not formally compared this in the current study, but a manuscript from our group related to this issue has recently been published (Preham et al CD4 +T cells alter the stromal microenvironment and repress medullary erythropoiesis in murine visceral leishmaniasis. Front. Immunol. 2018. Doi: 10.3389/fimmu.2018.02958).

4. Figure 6B: replace green square with green dots

This has been corrected

5. Figure 6E-F changes in granuloma cell populations over time

The F4/80 + DAPI images shown in Fig 1C-F were chosen to show representative granulomas of differing size and not to formally describe changes in cell populations over time. This type of analysis has been conducted previously e.g. Moore et al PlosOne 2012. 7(3):e34134 for B and T cells and Beattie et al PLoS Path 2010. 6(3):e1000805.   A further manuscript currently in preparation from the CRACKIT project (Hamp et al) will also contain further information on hepatic cellularity over time.  We have added an additional sentence to the discussion to highlight that transcriptomic changes as assessed here may be due to changes in cell composition or altered transcription levels within a given cell population.

6. Figure 10: comment on differential regulation of checkpoints in spleen and liver

We have included some additional discussion on these data and the difficulties in extrapolating to the effects of therapeutic blockade.

Tissue-specific immune responses are always of interest in visceral leishmaniasis. The authors address this question by looking at transcriptomes of the spleen, liver and blood of infected mice and comparing them with infected hamsters and human VL patients. I feel that their approach is reasonable and the data shown here are very useful to the research field. It will be great if the authors can strengthen the conclusion by discussing/interpreting the data more in depth.

Major comments

The authors described in the Discussion that BALB/c were used in this study ‘due to their extensive use as an immunological model of infection…’. However, the authors’ group also use C57BL/6 mice for immunological studies of experimental VL. Is there any stronger reason to select BALB/c over C57BL/6?In fact, pathogenesis/immune responses during L. donovani infection may be different between those two inbred mice. For example, in the current study using BALB/c mice, the splenic parasite burden and spleen weight were almost parallel. It means the parasite number per host nuclei are relative stable over the course of infection. In contrast, C57BL/6 mice show higher increase in LDU over spleen weight, according to Ato M et al., Nature Immunology, 2002. Also, splenomegaly by L. donovani infection seems worse in BALB/c mice (~8 times enlargement at D35) than C57BL/6 mice (only 1.7 time at D28). These differences should be welcomed because they mean that comparison of the inbred mouse strains is also a viable approach besides comparison with other species.

Although the title says tissue and host species-‘specific’ changes, I feel that the article more focuses on commonality among tissues and species. For example, the last part of this work focuses on the 26 genes commonly up-regulated in the spleen, liver and blood; what is the significance of identifying such genes? The authors know that at the time point kinetics of parasite proliferation in the spleen and liver are different, so what kind of phenomena should be associated with the commonly up-regulated genes? I believe that more careful watch into differently regulated genes, across tissues and across species, is indispensable. Identification of IDO-1 as a mouse-low, hamster-high and human-high factor is really intriguing and the authors are encouraged to take further analyses to identify such ‘protection-related factor candidates’ because it is supposed to be their objective.

Related to the previous comment, representation of each tissue/species is a little vague. My understanding is; mouse spleen, vulnerable at D42; mouse liver, protective at D42; hamster spleen, vulnerable at D42; hamster liver, vulnerable at D42; human blood, reflecting systemic vulnerable condition. Is it correct and if yes what does mouse blood represent? Although I understand the authors know that it is not that simple, data analyses and discussions with clear and consistent hypotheses/assumptions are helpful to the readers.

Page 4: Discussions on IL-10Ra and Trim24 are difficult to understand. ‘Upstream regulators that were predicted to be down regulated, despite not being DE themselves’, do you mean that they are upstream regulators for down-regulated genes? Also, what does this actually mean in the setting of murine VL? Can you also discuss the fact that IL-10 is up-regulated in the spleen at all time points and what it means in connection with IL-10Ra?

Not sure what Figure 3E means, what two kinds of line means there. No description in the legend. Besides the letters are too small to see.

Minor comments:

Many abbreviations are not spelled out at their first appearance; for example, DE and GSEA are explained in full at Results section while already used in Methods. Check other abbreviations too, some are not being spelled out at all.

Supplementary table for Figure 1A and 1B Raw data says the spleen weight of naïve mice at D36 is 412 mg. Please correct.

Page 9: The first sentence ‘Kong et al….’ please add the reference number here too.

I have read this submission. I believe that I have an appropriate level of expertise to confirm that it is of an acceptable scientific standard.

Major

1.  comment on strain choice and specific differences in response to infection

The reviewer is correct that we and others also use B6 mice extensively in studies of murine VL.  We have expanded in the Discussion on the rationale for our choice of BALB/c mice in this study.  However, we are currently preparing another manuscript from the CRACKIT program, which will include transcriptomic data from infected B6 as well as B6. Rag-/- mice.

2. why look at a common 26 gene signature? Differences across species should be pursued.

I think in this manuscript we have done both.  The value for a common gene signature might be in the identification of common underlying mechanisms or the identification of biomarkers that might be useful for clinical monitoring of disease progression.  We have expanded the discussion to address this.

3. clarify use of different time points in different species and relation to disease pattern

Unfortunately, we are limited by the availability of data.  Researchers working with hamsters have only studied spleen (and as their study was not published at the time we began our own investigations, we could not “time match” even if that was desirable). Researchers investigating human disease have to date only published transcriptomic data from peripheral blood. We hope our current manuscript will highlight the need for better comparative studies, but it will always be difficult to match to disease profile as these differ widely between models and indeed humans.

4. clarify discussion on upstream negative regulators

We have amended the text to explain how IPA determines upstream regulators

5. improve clarity and description of Figure 3E

We have included new text to clarify this figure.

Minor

1. check abbreviation call outs

These have been corrected

2. error in Figure 1A/B RAW data table

This has been corrected

3. page 9: add reference to Kong et al

This has been corrected