Tissue and host species-specific transcriptional changes in models of experimental visceral leishmaniasis [version 1; peer review: 3 approved, 1 approved with reservations]

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, hostpathogen 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 interferonregulated genes whilst confirming a greater activation of type 2 Open Peer Review


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 mantained 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 . 2x10 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 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. Comparison of the gene expression data between infected and naïve samples was performed using the Benjamini and Hochberg false discovery rate (FDR) correction 44 . 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 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. 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 . Venn diagrams were generated using Venny 2.1 47 . Hamster transcriptomic datasets were from 41. Human transcriptomic data sets were from 39.

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 2x10 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). To account for minor variations between groups of naïve mice used in the two experiments, differentially expressed (DE) genes were identified by reference to time-and experimentmatched 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 gene set enrichment analysis (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, hypoxia, myogenesis, allograft rejection, IL-6 -JAK/STAT signalling, angiogenesis, coagulation, G2M checkpoint, inflammatory response, oestrogen response, epithelialmesenchymal transition, E2F targets, KRAS signalling and complement. 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), despite not being DE themselves. For example, at day 36 p.i., activation of IFNγ was predicted by the regulation of 301 target genes, consistent with an observed Log2FC change in expression of 3.474. IFNγ was 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. The 41 genes that predict inactivation of TRIM24 during L. donovani infection are shown in Figure 3E.   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.
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   (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 ( Figure 5C). Amongst common down-regulated DE 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.
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).
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 analysis, 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. Of note, for many cytokines and chemokines, day 36 represented the peak of the response, 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 analysis 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 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, 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) negating quantitative comparison in expression levels, variables such as tissue parasite load and time of infection are difficult to standardise (particularly for human disease), and with the exception of our study, data is restricted to a single tissue site. 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 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, and clearly indicate the dominant nature of interferon signalling during murine VL. The clinical value of this signature is however questionable, 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 12 or differing analytical methods. It would be informative, given the above discussion, to compare human whole blood with hamster whole blood to ascertain whether the latter may provide a greater degree of similarity and to understand the extent of temporal and geographic diversity in the transcriptome of patients with VL.
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 CRACKIT 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).
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 . 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. Some of these checkpoint targets are equally represented in mouse and hamster spleen and human whole blood

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. The funders had no role in study design, data collection and analysis, decision to publish, or preparation of the manuscript.

Acknowledgements
The authors thank the staff at the Biological Services Facilities of University of York and LSHTM for their support with animal husbandry.

Is the work clearly and accurately presented and does it cite the current literature? Yes
Is the study design appropriate and is the work technically sound? 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? Yes
Competing Interests: No competing interests were disclosed.

I confirm that I have read this submission and believe that I have an appropriate level of expertise to confirm that it is of an acceptable scientific standard.
Author Response 17 Dec 2018 Paul Kaye, University of York, UK, York, UK

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.

issues resulting from comparison of RNA-Seq to microarray
We have expanded our discussion on this point as suggested by the referee

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

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

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

Is the work clearly and accurately presented and does it cite the current literature? Yes
Is the study design appropriate and is the work technically sound? 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
Competing Interests: No competing interests were disclosed.
We confirm that we have read this submission and 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.

Author Response 17 Dec 2018
Paul Kaye, University of York, UK, York, UK Methods

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

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