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.


Abstract
: Human visceral leishmaniasis, caused by infection with Background or is a potentially fatal disease affecting Leishmania donovani L. infantum, 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 infection in mice has been Leishmania donovani extensively characterized, transcriptomic analysis capturing the tissue-specific evolution of disease has yet to be reported. : We provide an analysis of the transcriptome of spleen, liver and Methods peripheral blood of BALB/c mice infected with . Where possible, we L. donovani 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 Our analysis confirms marked tissue-specific alterations in the Results: 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.
Our transcriptomic data, coupled with histopathologic analysis of Conclusions:

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][16][17][18][19][20][21][22][23][24][25][26] , and on the response to drugs, vaccines and other immunotherapeutic interventions [27][28][29][30][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][36][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][41][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.

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

Amendments from Version 1
We thank the reviewers for their constructive and supportive comments, which we have fully addressed in the specific responses. In summary, we have i) improved the clarity of some of the Figures, ii) added more technical details into some sections of the Materials and Methods and to the Results sections where bioinformatic analysis is discussed, iii) added further discussion on possible batch effects, the choice of animal model for our transcriptomic analysis and the potential pitfalls of comparing data derived from multiple platforms, and iv) corrected a transcribing error in the RAW data file. Minor editorial corrections have also been made to improve readability. Finally, as this study is part of a larger program of work (the CRACK IT Virtual Infectious Disease Research project), we have also now included for the readers benefit, a brief description of this project and an indication of related manuscripts that are currently in preparation for publication.

REVISED
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 2 O. 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 2 O 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 2 FC = 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 2 FC -3.52 to +5.82; d17: 2826 probes, log 2 FC -4.01 to +6.37; d21: 5124 probes, Log 2 FC -4.54 to +7.59; d36: 5051 probes, Log 2 FC -6.45 to +8.36; d42: 4946 probes, Log 2 FC -6.30 to +7.13).
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   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. 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).
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.  (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).
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.
Of the top 25 down-regulated genes, 10 encoded amylase or amylase precursor proteins and 10 have serine peptidase activity (Table S3).
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  Table S3.  process of granulomatous inflammation has been wellcharacterised 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).
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.
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 .
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).
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).
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).
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).

Discussion
The immune response to L. donovani infection has been extensively analysed using a wide array of technologies, including   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 bromodomaincontaining 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.

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 Figure 3A.   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. 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.

If applicable, is the statistical analysis and its interpretation appropriate? Yes
Are all the source data underlying the results available to ensure full reproducibility? Yes

Are the conclusions drawn adequately supported by the results? Yes
No competing interests were disclosed. Competing Interests:

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.
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.

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.

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).

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.
No competing interests were disclosed. 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 1.

4.
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.

Is the work clearly and accurately presented and does it cite the current literature?
Yes 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".

Is the study design appropriate and is the work technically sound? Yes
Are sufficient details of methods and analysis provided to allow replication by others? Yes

If applicable, is the statistical analysis and its interpretation appropriate? Yes
Are all the source data underlying the results available to ensure full reproducibility? Yes

Are the conclusions drawn adequately supported by the results? Partly
No competing interests were disclosed.

Competing Interests:
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.

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.

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.

Fig1A/B Raw data on spleen size in naïve controls at d36 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

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). Figure 6B: replace green square with green dots This has been corrected 5. Figure  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. 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.

6.
No competing interests were disclosed.

Yasuyuki Goto
Laboratory of Molecular Immunology, Graduate School of Agricultural and Life Sciences, The University of Tokyo, Tokyo, Japan 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. + Figure 3E We have included new text to clarify this figure.

check abbreviation call outs
These have been corrected 2. error in Figure 1A