Splenic Transcriptional Responses in Severe Visceral Leishmaniasis: Impaired Leukocyte Chemotaxis and Cell Cycle Arrest

Structural changes in the spleen have been reported in several infectious diseases. In visceral leishmaniasis (VL), a severe parasitic disease caused by Leishmania spp., the loss of white pulp accompanies a severe clinical presentation. Hamster model reproduces aspects of human VL progression. In the early stages, a transcriptomic signature of leukocyte recruitment was associated with white pulp hyperplasia. Subsequently, impaired leukocyte chemotaxis with loss of T lymphocytes in the periarteriolar lymphoid sheath occurred. This differential gene expression was subsequently corroborated by transcriptomic profiling of spleens in severe human VL. At the latest stage, spleen disorganization was associated with increasing clinical signs of VL. White pulp disruption was accompanied by decreased DLK1 expression. The expression of CXCL13, CCR5, CCL19, CCR6, CCR7 and LTA decreased, likely regulated by CDKN2A overexpression. Our findings enlighten a pathway implying cell cycle arrest and decreased gene expression involved in spleen organization.


INTRODUCTION
The spleen plays a central role in the pathogenesis of visceral leishmaniasis (VL), a severe parasitic disease caused by Leishmania infantum or Leishmania donovani (1)(2)(3)(4). In VL, infection in the spleen is persistent. This organ presents important sequential morphological changes that parallel the general clinical state of disease progression (5)(6)(7). The lymphoid follicles located in splenic white pulp (WP) respond to Leishmania infection by forming the germinal center and increasing B cell differentiation. This hyperplasia in the WP is accompanied by increasing numbers of Leishmania-infected macrophages, which leads to splenomegaly and increased parasite burden in later stages of disease (8). The enlargement of the spleen exacerbates its function of blood cell recycling, contributing to pancytopenia (9,10). As the disease progresses, cell populations become substantially altered in the white and red pulp (RP). Lymphoid follicles decrease in size and the boundaries between the germinal center mantle area and marginal zone become less evident, eventually disappearing (6,(11)(12)(13).
Disruption in splenic microarchitecture has been reported in a variety of diseases, such as leishmaniasis and malaria, as well as many viral diseases, including HIV and COVID-19 (14)(15)(16)(17). However, the mechanisms underlying these changes have yet to be elucidated. The splenic architecture is chemically maintained by a complex network of transcripts. In Leishmania-infected dogs with disorganized spleens, decreased CXCL13 expression has been associated with a loss of follicular dendritic cells and decreased numbers of B cells in lymphoid follicles (18). Concomitantly, the accumulation of plasma cells in the RP may result from interference in B cell differentiation pathways in the spleen (13). Moreover, the overexpression of CXCL12, BAFF and APRIL in splenic RP supports the hypothesis of anomalous homing and increased plasma cell life span (13). In the periarteriolar lymphoid sheath (PALS), T cell apoptosis driven by cellular exhaustion markers, such as CTLA-4, may lead to decreased numbers of CD4+ T lymphocytes, which has been associated with WP disorganization (18)(19)(20).
Hamsters are known to develop a progressive form of VL that can be fatal, recapitulating many aspects of the disease observed in susceptible dogs and humans (21,22). In addition to being associated with a mixed Th1 and Th2 response profile, susceptibility has also been linked to a failure to express the gene coding for nitric oxide synthase (iNOS), even in response to IFN-g, which implies the inability of macrophages to control replication of intracellular parasites (23,24). Although experimental infection in hamsters represents a suitable model for the study of lesions associated with severe forms of VL, the analysis of pathways involved in lesion development has been limited by a lack of compatible reagents. However, recent gene expression studies have substantially contributed to our understanding of mechanisms involved in immune response to VL (25,26).
The present study demonstrates that L. infantum infection leads to disorganization of the spleen compartments in hamsters with severe clinical presentations of VL. We employ large-scale gene expression analysis to identify transcripts potentially involved in WP disruption. We propose that the histological changes seen during the course of infection correlated with the transcriptional differential gene expression and tissue distribution of T lymphocyte marker and DLK1 protein. We then corroborate our findings from experiments performed in hamsters using splenic transcriptomic profiles of severe human forms of VL.

Ethics Statement
All experimental procedures were approved by the Institutional Review Board of the Goncalo Moniz Institute (IGM-FIOCRUZ, license nos. 004/2013 and 017/2015 for experiments performed in hamsters, and 3.491.092 for humans), and were carried out in accordance with Brazilian legislation on ethical animal experimentation. All patients were informed and consented to participate in the study.

Hamsters and Infection Procedures
Golden Syrian hamsters (Mesocricetus auratus) were obtained from a FIOCRUZ animal care facility. Male hamsters, 6-8 weeks old, weighing 114 ± 20.2 g, were kept in cages under a controlled physiological light and temperature conditions. Each hamster was individually identified by a subcutaneous chip that was recognized by compatible reader (Microchip Partners). Injections were performed by intraperitoneal route in 1mL of either saline solution (control group) or a parasite suspension (infected group) containing 10 7 promastigotes of L. infantum (strain MHOM/ BR2000/Merivaldo2). L. infantum virulence was maintained through passages in hamsters, and parasites were grown until reaching stationary phase in vitro in Schneider's complete medium (Schneider + 20% fetal bovine serum [FBS], Gibco, USA) in a B.O.D. incubator at 24°C. Hamsters were euthanized by anesthetic overdose (10mg cetamin + 1mg xylazine/mL) at 30, 60, 120 and 150 days post-injection (dpi). Animal spleens were collected, weighed, and fragmented into three pieces for histopathological analysis, parasite culturing and RNA extraction. For parasite detection, macerated spleen fragment was sown in biphasic culture medium (blood agar + complete Schneider medium) in a B.O.D. incubator at 24°C and examined for up to 4 weeks.

Clinical Signs, Serum Biochemistry and Hematology
Clinical signs of VL (presence of skin lesions, hair and weight loss, splenomegaly) were examined on a weekly basis. At the day of euthanasia, blood was collected by heart puncture. Blood was preserved in EDTA blood collection tubes and submitted to analysis of total red blood cell (RBC) counts and total and differential white blood cell (WBC) counts using an automated cell counter. For biochemical analysis, the blood was collected in serum collection tubes and centrifuged at 2,000 xg for 15 minutes at room temperature and serum was submitted to biochemical analysis of total serum protein, albumin and globulin fractions.

Human Spleens
Human spleen samples were obtained from the Natan Portella Tropical Disease Institute (Teresina, Piauı, Brazil) following therapeutic splenectomy in patients with severe VL, two of whom were also coinfected with HIV. Both coinfected VL-HIV patients had undetectable viral load. Splenic material from two males and one female, median age of 51 [range: 45 -55] years, were included in the study. Spleen fragments were submitted to histopathological analysis and RNA extraction. in hematoxylin & eosin. All tissue specimens were examined by two pathologists (WLCdS and LARF) on a group-blinded basis to estimate the intensity of inflammatory infiltrate, fibrosis and cell degeneration, as previously described (27). Architectural organization was assessed in the spleen samples according to Hermida et al. (3). In brief, well-organized or type 1 spleen was considered when WP microenvironments were easily distinguished in lymphoid follicles, the germinal center, PALS, mantle zone and marginal zone; slightly disorganized or type 2, corresponded to the loss of some boundaries between the WP microenvironments due atrophic or hyperplastic alterations; extensively disorganized, or type 3, was considered when boundaries were poorly delimited between WP and RP regions, and lymphoid follicle atrophy was observed (3,6). Morphometric evaluations were performed to estimate WP/RP ratios using ImagePro Plus (Media Cybernetics).

Immunohistochemistry and Morphometry
Spleen fragments were fixed in 10% formalin and embedded in paraffin; 4 mm-thick sections were obtained and mounted on silanized slides. The slides were then deparaffinized and rehydrated. For CD3 staining, the slides were immersed in a 10% ammonium hydroxide solution for 10 minutes to remove formaldehyde pigment, 3% hydrogen peroxide for 10 minutes to block endogenous peroxidase and in Tris-EDTA pH 9.0 for 30 minutes at 96°C to perform antigenic recovery. Slides were then cooled to room temperature for 20 minutes, and nonspecific staining was subsequently blocked by 5% bovine serum albumin for 20 minutes. For DLK1 staining, the slides were incubated with 20% ammonium hydroxide in 95% alcohol for 20 minutes, washed in running tap water, followed by washing in distilled water. Antigen recovery was performed by heating (115°C) in citrate buffer (pH 6.0) for 20 minutes. Endogenous peroxidase was blocked with 3% hydrogen peroxide in 3 incubations of 10 minutes, followed by washing with distilled water. Unspecific staining blocking was performed with horse serum (species of the secondary antibody) (Vector, 30022) for 15 minutes and washed in PBS, followed by permeabilization in PBS with 0.1% saponin for 5 minutes in 2 incubations. Afterwards, samples were incubated overnight with either anti-CD3 (1:300, ab16669, Abcam), rabbit anti-DLK1 (1:1000, ab21682, Abcam) or anti-CD20 (1:200, ab64088) primary antibodies, followed by incubation with the secondary antibody conjugated HRP (ab205718, Abcam) for 30 minutes for CD3 and CD20 staining, or anti-rabbit IgG polymer (Vector, 30026) for DLK1 staining. Staining was developed using a 0.02% 3,3-diaminobenzidine solution for 5 minutes, followed by nuclei counterstaining in Harris hematoxylin (Sigma). Spleen sections for CD3 staining were scanned at 200x magnification using an Olympus VS110 virtual slide scanning system (Olympus America Inc, United States). CD3 measurements were estimated in the five largest WP regions in seven cases per group, without overlapping areas. The boundaries of WP and PALS were established according to previously described morphological criteria (28,29), and confirmed by staining for CD20 to visualize lymphoid follicles in WP. Images of the three hotspot areas of DLK1 staining were obtained at 400x magnification using Olympus BX53 camera in Image-Pro Plus software (V. 4.5, Media Cybernetics) and cell density was estimated. The same spleen specimens submitted to transcriptome analysis were used for DLK1 immunohistochemistry. Unfortunately, the specimens of 30 dpi timepoint, had a strong background that impaired the morphometric estimation. All measurements were morphometrically estimated using ImageJ software version 1.52 (National Institutes of Health, United States) following the delimitation of each respective region of interest. Additionally, the total area of CD3 staining was estimated in each compartment. DLK1 staining morphometry was estimated by two observers and the mean between the measures was represented. The means of cell density were used for statistical analysis.

RNA Isolations
Spleen fragments collected from hamster dissection and biopsied humans were stored at -80°C and dissociated in Trizol (Invitrogen). Total RNA isolation was performed using the aqueous phase of Trizol mixture followed by purification using RNAeasy extraction kit (Qiagen, Inc.) according to the manufacturer's protocol, including an extra step of DNA lysis through the addition of DNase. Total RNA purity was determined by NanoDrop spectrophotometer (ThermoFischer) and the concentration of the material was measured in Qubit fluorimeter (Invitrogen). Total RNA integrity was evaluated on an automatized electrophoresis Agilent Bioanalyzer 2100 (Agilent) with RNA 6000 Nano kit (Agilent). Samples were considered of high quality when RNA integrity number (RIN)>7.

Hamster Spleen Transcriptomics
Large-scale expression data was obtained to identify the transcriptional response, including differentially expressed genes, between infected and non-infected animals in different time-points (30, 60, 120 and 150 days). Additionally, L. infantum RNA counts were also assessed. Once a total RNA was obtained that met the purity, concentration and integrity criteria, a contaminant depletion step (globins and ribosomal RNAs) was applied and the resulting sample was converted into cDNA and subjected to NextSeq platform (Illumina, Inc), according to the manufacturer's protocol. An average of 30,000,000 readings per sample was expected. The sequencing was performed using pairedend strategy with fragment library of size 2×75 bases using stranded protocol. Transcriptomic profiles from the experimental VL samples were obtained using three biological replicas for each group (control and infected) at each of the four timepoints studied. The sample selection for transcriptomic analysis was based on the disruption of the splenic WP throughout the followup of the infection. The quality of reads was assessed by FastQC (v0.10.1) in each sample, followed by summarization of the reports using MultiQC (30). Golden hamster (Genome ID 11998) or Leishmania infantum (Genome ID 249) reference genomes were obtained from the NCBI Genome database and then aligned using STAR (31). The number of reads per transcript were calculated by FeatureCounts. Gene annotation was performed and differential expression (DE) was determined using data from the GEO repository (NCBI). Differential expression analysis was conducted using edgeR (32). Wald's test was performed using DESeq2 (25). Functional analysis was carried out using Ingenuity Pathway Analysis (IPA) software (Qiagen, Inc.) considering a cut-off log fold change de Melo et al.
Transcriptional Responses In Visceral Leishmaniasis (FC) value of 2, false discovery rate (FDR) <0.05 and absolute z score of ≥2 (26). Unsupervised analysis was employed at each timepoint for the evaluation of canonical pathways and enrichment of diseases and functions (Tables S1, S2). Enrichment analysis of differentially expressed genes was performed overlapping networks for each timepoint. To investigate transcripts associated with spleen morphology and WP organization, a supervised analysis was performed with the addition of a panel of transcripts, followed by an unsupervised network expansion to evaluate interaction pathways related to splenic organization.

Human Spleen Transcriptomics
Gene expression was assayed using nCounter platform (NanoString Technologies, Seattle, Washington), based on direct molecular bar codes coupled to target RNA transcripts using digital detection at the Genomics Core Leuven (VIB/KULeuven -Belgium). Two probes were employed, one that captures the mRNA of interest for the complementary sequence, while the other is connected to a fluorescent barcode that identifies the target by specific hybridization. The probes are mixed with the genetic material and transcripts of interest are identified by fluorescent barcodes. Samples from healthy spleen tissue were reanalyzed and obtained from the GTEx project (33), which performed whole-transcriptome profiling of multiple human tissues, including spleen. Three samples were included for this analysis, and raw reads from whole-transcriptome profiling were obtained from the SRA selector at http://trace.ncbi.nlm.nih.gov/Traces/study/?acc= phs000424. These samples were chosen based on their age and gender profile similarity with our infected tissue human donors. Gene annotation was performed using the biomaRt library for R (34). Mapping of the raw reads against the human genome reference (GRCh38) was performed using STAR (31) at default parameters, and HT-seq count was used to assign and count mapped reads to annotated genomic features (35) using GENCODE v. 25 annotations. Nanostring expression was quantified by counting the amount of mRNA identified in a single reaction. The present study constructed a panel consisting of 9 genes for human analysis (Table S3). Raw data were preprocessed using nSolver 2.0 software (NanoString Technologies).

Statistical Analysis
The statistical analysis of transcriptomic data was performed as described above. Absolute numbers, means, medians and percentages are indicated in the text, tables and graphs. Comparisons between control and infected groups were performed using the t test or Mann-Whitney test as indicated.
To perform comparisons between control and infected groups at each time point, ANOVA was used in combination with multiple comparisons by the Kruskal-Wallis test. Statistical significance was considered when p<0.05.

Clinical Evaluation
The general characteristics of the animals are presented in Table 1. All infected animals presented positive spleen culture throughout all time points, and RNA copies aligned to Leishmania infantum presented an increasing trend over time of infection. Clinical signs of disease were present after 120 dpi in two animals (one had a crust on the snout, and another presented an ulcer in the oral region) and more intense in two animals at 150 dpi (one had cachexia, irregular breathing, dehydration, an ulcer in the oral region and pleural effusion; another was emaciated and presented an ulcer in the oral region). Statistically significant differences in clinical presentation were observed between the groups at 150 dpi ( Table 1).
While weight was consistently lower in infected animals compared to controls, significant differences were only seen at 30 and 150 dpi. Despite an absence of statistical significance, differences were detected at 150 dpi in hematological parameters, with decreasing trends observed in average red blood cells (RBC) counts, hemoglobin, hematocrit, RBC mean corpuscular volume (MCV) and white blood cell (WBC) counts in the infected animals. Upon necropsy, splenomegaly was observed at 120 and 150 dpi ( Table 2).

Histological Evaluation of Spleen
Chronic perisplenitis was observed in 4/7 infected animals at 150 dpi ( Table 2). Smaller proportion of WP were seen in the L. infantum-infected animals compared to controls after 60 dpi, with statistically significant differences observed between the groups at 120 dpi. At 150 dpi, despite a persistent trend towards decreased WP size in the L. infantum-infected group, controls also presented less WP, possibly due to aging ( Table 2). The disruption of WP architecture was observed in 10 L. infantum-infected hamsters, considered slight (Type 2) in 5/7 (71.4%) animals at 120 dpi and in 3/7 (42.8%) animals at 150 dpi, and intense (type 3) in 2/7 (28.6%) hamsters at 150 dpi (Mann-Whitney test, p=0.03) ( Table 2 and Figure 1). WP disruption was associated with germinal center atrophy in 3/7 (43%) infected hamsters at 150 dpi ( Table 2).

Transcriptomic Splenic Profiles
Based on the classification of WP organization, spleen tissues from uninfected (n=12) and infected (n=12) hamsters were submitted to RNA-seq analysis. Overall, the transcriptomic profile revealed an increase in differentially expressed (DE) transcripts during the course of infection ( Figure S1A): 1040 (30 dpi); 1813 (60 dpi); 5874 (120 dpi) and 6512 (150 dpi). Principal component analysis pictures differential distributions in accordance with infection status and duration ( Figure S1B). The distribution of groups along PC2 reflects some degree of splenic disorganization, while organized spleens are distributed along PC1 ( Figure S1B). To investigate the functional effects of smaller groups of genes, an adjusted false discovery ratio (FDR) <0.05 and fold change (FC) of 2 were used to generate consistently expressed interaction networks at each time point. Under these parameters, only two transcripts were included in the analysis at 30 dpi, reflecting a lack of significant alterations in the histology of the spleen and the absence of clinical presentation at 30 dpi. However, at all remaining time points, higher numbers of transcripts were included making it feasible to build interaction networks. At 60 dpi, a set of overexpressed genes were predicted as regulator effectors for the functions of leukocyte trafficking and inflammatory process (Figure 2A). At this timepoint, genes related to granulocyte recruitment were also found to be overexpressed. At 120 dpi, the granulocyte recruitment remained activated, with a greater number of genes included ( Figure 2B). At this time point the quantity of ion metal became altered, primarily stimulated by SPP1, CCL5, IFNg, LEP and CXCL10. This effect was enriched to the last time point (150 dpi), revealing an increased number of associated transcripts ( Figure 2C). Systemic alterations as in hemorrhagic disease and behavior were predicted in response to the set of transcripts expressed at this stage of the infection. Part of the differential expression of transcripts involved in the functional analysis was shared between time points 60 and 120 dpi ( Figure 2D) and between 120 and 150 dpi ( Figure 2E), reflecting a kinetics in the molecular signature.

Transcriptomic Kinetics Processes
Given the importance of analyzing biological processes along the course of infection, we overlapped transcripts identified as differentially expressed onto enriched networks to elucidate key processes (Figure 3). At 30 dpi, as only two differentially expressed transcripts were identified after the functional analysis (FDR<0.05, FC of 2), only Delta Like Non-Canonical Notch Ligand 1 -DLK1 was observed to be differentially expressed in common with the remaining timepoints. At 30 dpi, DLK1 was observed to be upregulated ( Figure 3A, 30 dpi). This expression pattern was reversed at 150 dpi (the latest stage of infection), when the downregulation of DLK1 occurred ( Figure 3A, 150 dpi). We performed an antibody-based validation of this transcript of DLK1 in the splenic tissue at 60, 120 and 150 dpi ( Figure 3B). Cells expressing DLK1 were distributed around perivascular area. The density of DLK1expressing cells decreased at each timepoint (p=0.0002), correlating with the decreasing expression of DLK1 mRNA found in the transcriptomic analysis ( Figure 3C) of the spleen of Leishmania-infected hamsters. With regard to gene expression analysis by overlapping along infection, transcripts deemed relevant to VL, such as ARG1 and IDO1, exhibited high expression at 60, 120 and 150 dpi ( Figure 3D).

T lymphocytes and Spleen Disorganization
Splenic T lymphocytes identified by CD3 staining were found on decreasing trend in area percentages of PALS in infected hamsters from 60 dpi (Figures 4A, B). Strikingly, when we performed an unsupervised enrichment of biological processes, chemotaxis of T lymphocytes emerged with the highest score in the dataset and consistently presented a pattern of differential expression over the course of infection ( Figure 4C). The upregulation of the T lymphocyte chemotaxis pathway observed at 60 dpi decreased to the next timepoint (120 dpi), concomitant with decreases in PALS and WP areas. Impaired T lymphocyte chemotaxis persisted until 150 dpi and in association with statistically significant lower PALS area, indicating a potential transcriptomic mechanism linked to profound WP disorganization at later stages of disease. Genes associated with T lymphocyte chemotaxis included SPP1, IL-21, CXCL9, IFN-g, CXCL11, CCL5, CCL24 and CXCL10, which exhibited upregulation, as well as CCL21 and IL12B which appeared downregulated ( Figure 4D). Except for CCL17 and CCL21 that were not included in the nCounter target panel, all the transcripts involved in T lymphocyte chemotaxis were validated in the splenic transcriptome in human VL, of which included two organized (type 1) and one disorganized spleen (Type 3) ( Figure 4D). Relevant correlations between tissue expression of CD3 and transcripts of the genes in the T lymphocyte chemotaxis pathway were observed, indicating a genetic regulation of splenic organization ( Figure 4E).

Cell Cycle Arrest and Spleen Disorganization
In light of concurrent findings regarding impaired leukocyte recruitment and the histological disruption of spleen microenvironments, we then investigated expression patterns among the following set of genes known to participate in spleen compartment organization: CXCL13, CXCR5, LTA, CCR7, CCR6, CCL19, CCL21, CXCL12 and MKI67. To evaluate transcripts associated with spleen morphology and WP organization, we employed unsupervised network expansion to evaluate interaction pathways related to the function of splenic organization. After enrichment with 20 unknown connecting transcripts, a network associated with the regulation of spleen disorganization was constructed. While Cyclin dependent kinase inhibitor 2 A (CDKN2A) was found to be upregulated at 150 dpi ( Figure 5A), CXCL13, CCR5, CCL19, CCR6, CCR7 and LTA were downregulated, suggesting an inhibitory role played by CDKN2A. Interestingly, in our unsupervised analysis of canonical pathways, the consistent upregulation of CDKN2A indicated its involvement in the down-modulation of cell cycle regulation ( Figure 5B).

Infection by L. infantum Induced Morphological Changes Leading to Spleen Disorganization in Hamsters
Several experimental animal models of infection have been used to study changes brought about by VL. In mice and hamsters, Veress et al. (11) observed the hyperplasia of lymphoid follicles during the progression of disease (11). In BALB/c mice, substantial tissue remodeling and changes in cell populations were found following infection by L. infantum or L. donovani (22,27,36,(38)(39)(40). However, the marked changes of the WP observed by other authors at early stages of infection did not progress to WP disruption at later stages of infection in mice (22,27,36,(38)(39)(40). In contrast, hyperplasia followed by atrophy and white pulp disruption was observed at later timepoints in the present experimental model involving hamsters. Veress et al. (22) observed complete WP disruption at around 120 days of infection, accompanied by extensive RP replacement by amyloid deposition (22). Herein, intense WP disruption was observed  only at 150 dpi in the absence of amyloidosis. These differences may be associated with the Leishmania strain and inoculum size used in experimentation (38). Despite the sample size of seven animals per group, and the fact that hamster model is susceptible to severe forms of VL, not all the infected animals developed complete disruption of the WP during the time course of this study. In fact, in canine and human VL not all the individuals develop spleen disorganization, as reported by Veress et al. (22) and Santana et al. (6). Nevertheless, the L. infantum infection successfully induced splenic disorganization in 5 out of 7 animals at 120 dpi, and 3 out of 7 at 150 dpi. Further, 2 out of 7 animals had extensive disruption of the WP at 150 dpi. Our data indicate more exacerbated clinical VL manifestations at later stages of disease, as evidenced by weight loss, splenomegaly and a greater frequency and intensity of clinical signs.

Transcript Expression Kinetics Uncovers DLK1 Deregulation
The clinical and histological changes identified herein were found to be correlated with different transcriptomic signatures in accordance with the progression of infection. A normal resting spleen presents with balanced expression of genes involved in the   (8,12,36). Along the course of the present experimental infection model, we observed a progressive decrease in DLK1 expression, with WP disruption taking place by 150 dpi. Immunohistochemical staining for DLK1 corroborated the transcript expression and spleen disorganization. We chose an antibody-based validation for this finding, to reflect its biological relevance in the protein synthesis. As we show herein, cells expressing DLK1 decreases in density in the spleen over the course of infection. The DLK1 gene is a member of a family of growth factors that regulate cell differentiation, and is highly expressed by CD34+ hematopoietic cells (41). DLK1 has also been implicated in the differentiation of several cell types regulated by the Notch gene family (42). The absence of DLK1 expression in DLK1 -/mice was shown to lead to decreased frequency of B cells in the lymphoid follicles, an increased presence of B cells in the MZ and greater IgG1 production (37). We therefore speculate that losses in follicular B cell populations, the accumulation of IgG-producing plasma cells in the RP and associated hypergammaglobulinemia previously reported by other authors may in fact be associated with alterations in DLK1 expression as observed herein during the progression of infection (6,13). Downregulation of DLK1 was reported in mice infected with L.
donovani at 21 and 42 days of infection (26). Although no correlation with spleen remodeling was investigated, it was shown that morphological changes in the spleen of mice with VL take place as early as 14 dpi. However, the spleen disorganization seems not to progress at later time points of the infection (8,12,27,39,40). In sum, correlations between changes in DLK1 expression and the remodeling of splenic microenvironments seem to suggest the participation of this transcript in the regulation of B cell differentiation, likely leading to WP disruption and possibly even the accumulation of plasma cells in the RP.

Splenic Transcriptomic Signatures Are Associated With Biological Functions of Leukocyte Trafficking and Granulocyte Recruitment
Transcriptomic processes involved in inflammatory response, leukocyte migration, granulocyte recruitment and lymphocyte chemotaxis were identified concomitantly with the detection of granulomas at 60 dpi. The establishment of infection at this stage could occur in response to the high expression of IDO1 and ARG1. IDO1 is induced by IFNg and has been associated with the polarization of M2 macrophages in hamsters with VL (25). In Leishmania infection, IDO1 also shows a function of suppressing innate immunity through regulatory stimulation of dendritic cells, and adaptive through the generation of regulatory T lymphocytes (43,44). ARG1 expression is induced in macrophages by Leishmania, which competes with iNOS for arginine and produces nutrients that favor parasite survival, such as urea, as well as decreasing the production of parasitotoxic nitric oxide (45). Altogether, these observations provide evidence of an inflammatory process at 60 dpi, which likely induced the T lymphocyte chemotaxis observed in the transcriptome at this stage. T lymphocyte chemotaxis is notably high between 30 and 60 dpi, and granulocytes recruitment is increased, associated with differential expression of the genes presented in the Figures 2, 4.
Previous data from our group demonstrated that normal RP cellularity becomes replaced by high numbers of plasma cells in association with the elevated expression of BAFF, APRIL and CXCL12 (3,13). In fact, it is possible that granulocyte recruitment may be related to the presence of eosinophils in the splenic compartment, which could supply factors promoting the survival and maintenance of plasma cells at later times of infection (46). A significant downmodulation in lymphocyte chemotaxis associated with the decreased expression of IL12b and CCL21 was identified at 120 dpi, which correlated with the observation of moderate disorganization of splenic lymphoid tissue. Other authors previously reported a highly inflammatory environment in the spleen of hamsters infected by L. donovani (25). IL12b expression is associated with the polarization of a Th1 response against pathogens in the spleen, whereas CCL21 expression is important for the retention of T lymphocytes in the periarteriolar zone of WP (29,47). Alterations in the expression of these transcripts were correlated with decreasing tissue expression of CD3 in the PALS area, supporting a genetic regulation of splenic placement of T lymphocytes in the progression of the disease. Besides providing a validation of the alteration on splenic T lymphocytes found in the transcriptome at the cellular level, we further sought to provide insights on translation of the hamster model to the human disease. Using the nCounter method for deep-genome and multiplexed profiling, the same molecular signature was found in human VL. At 150 dpi we found profound disorganization of the spleen and decreased expression of transcripts involved in the organization of splenic microenvironments. The reduction of PALS in splenic tissue has been associated with decreased expression of CCL19 and CCR7, the ligand of both CCL19 and CCL21 (40,48). We also detected a progressive reduction in the expression of LTa, CXCL13 and its receptor, CXCR5, along the course of infection. LTa is produced by stromal cells and is involved in WP formation during embryogenesis (29). In mice, the development of VL was marked by decreased number of stromal cells and impaired migration of dendritic cells within spleen compartments, both associated with low expression of CCL19 and CCL21 and attributed to a defective function of CCR7 (48). In canine VL, despite the absence of alterations in CCL19 and CCL21 expression, splenic disorganization was found to be associated with a decreased frequency of CD3 T lymphocytes (18), T lymphocyte apoptosis (19) and T cell exhaustion (20). Moreover, a reduction in B cell density in lymphoid follicles was accompanied by a decreased CXCL13 expression in dogs with VL and disorganized spleens (18).

Transcriptomic Signature of Spleen Disorganization Is Associated With Cell Cycle Arrest
At the final stage of infection (150 dpi), the observed differentially expressed transcripts formed a network of gene regulation associated with spleen disorganization. When an unsupervised network expansion involving these transcripts was performed, a central link was identified with CDKN2A. Throughout the kinetics of gene expression shown herein, CDKN2A expression was observed to change from reduced (30 dpi) to an increased state at all timepoints after 60 dpi, coinciding with the development of an immunosuppressive environment in the spleen that resulted in poor immune response observed against L. infantum, and possibly other pathogens (49). Correspondingly, the expression of cytokines, chemokines and receptors associated with disorganization (CXCL13, CCR5, CCL19, CCR6, CCR7 and LTA) predicted regulatory function by CDKN2A and was observed to significantly decline after 60 dpi. Immunohistochemistry staining corroborated the impairment of T lymphocyte chemotaxis between 120 and 150 dpi, leading to a reduction in PALS area, which may be potentially associated with the downregulation of CDKN2Arelated genes as shown at 150 dpi. These changes are reflected on a tissue level, as evidenced by the fading of MZ, reduced PALS area and lymphoid follicle atrophy.
CDKN2A is known to transcribe information used in the production of proteins, such as p16 (INK4A) and p14 (ARF), both involved in cell cycle inhibition in the G1 phase (50). This transcripts also plays a role in cell senescence and is involved in immune inflammatory events leading to type I diabetes and atherosclerosis, in addition to being associated with inflammatory bowel disease (51,52). Taken together, the associations reported herein suggest a regulatory role for CDKN2A with respect to cytokine-producing cells. CDKN2A was previously found upregulated in both hamster and mouse experimental models of L. donovani infection at earlier time points in comparison to this study (25,26). Importantly, this transcript has been extensively investigated as a biomarker in some types of cancer (50), and may also constitute an early marker of severe outcomes in VL.
The present work explores transcriptional pathways involved in leukocyte chemotaxis, cell cycle arrest and possible alterations in B cell differentiation, in association with histological and clinical changes during the course of VL. The differential gene expression correlated with the phenotypical changes in the T lymphocytes in the spleen, and a similar pattern of the expression of transcripts involved in the spleen remodeling was found in severe cases of human VL. The transcriptome of human VL is, however, limited due to the lack of paired controls. The acquisition of healthy splenic tissue is difficult, since this organ is not usually submitted to biopsy. In some cases it is obtained from splenectomy due to trauma. We made use of the GTEx database (33), which is well established and used by many others. Differential expression of DLK1 and CDKN2A emerging with potential regulatory roles in spleen disorganization is novel in the field and requires further investigations to confirm its relevance.
Our results serve to affirm, on a transcriptomic and histological level, the central role played by the spleen in the establishment of infection and disease progression and stand in agreement with previous phenomenological observations associating a severe profile of VL with disorganization of the spleen.

DATA AVAILABILITY STATEMENT
The datasets presented in this study can be found in online repositories. The names of the repository/repositories and accession number(s) can be found below: https://www.ncbi. nlm.nih.gov/, PRJNA695421.

ETHICS STATEMENT
The studies involving human participants were reviewed and approved by Instituto Goncalo Moniz IGM-FIOCRUZ. The patients/participants provided their written informed consent to participate in this study. The animal study was reviewed and approved by Instituto Goncalo Moniz.

AUTHOR CONTRIBUTIONS
CM, MH and WS designed the animal experiments. CM, MH, JF, BM, and RB performed the experiments. RK, FG, PR, and GF designed the transcriptomic study and RK, FG, and CM conducted the analysis. LF, CC, and WS obtained the human samples and conducted the histopathological analysis. CM and WS wrote the manuscript. All authors contributed to the article and approved the submitted version.

FUNDING
The authors acknowledge the financial support provided by the Conselho Nacional de Desenvolvimento Cientıfco e Tecnoloǵico -CNPq (grant no. 400905/2013-2) and the research productivity scholarship provided to WS. We are also thankful for financial support from the Coordenacão de Aperfeicoamento de Pessoal de Nıvel Superior -CAPES (grant no. 2072_2013) and the Fundacão de Amparo à Pesquisa do Estado da Bahia -FAPESB (no. PET0053_2013).