M199, FBS, Molecular grade water, methanol, chloroform, ribitol (adonitol), methoxyamine HCl, pyridine and primers from Sigma-Aldrich, St. Louis, MO, USA. BSTFA was obtained from Supelco Sigma, Bellefonte, Pennsylvania, United States. SYBR™ green PCR master mix was obtained from Applied Biosystems Foster City California, United States. Ketamax 50 was obtained from Troikaa Pharmaceuticals Ahmedabad, Gujrat, India. DNeasy blood and tissue kit was provided by Qiagen Hilden, Germany. A strain of L. donovani, BHU1260, was derived from the splenic aspirate of a VL patient at the Kala Azar Medical Centre of the Institute of Medical Sciences, Banaras Hindu University, Varanasi, India. L. donovani strain was grown in M199 medium (pH 7.4), which was supplemented with 10% fetal bovine serum (FBS) (Gibco, Carlsbad, CA, USA) and maintained at 23°C.

All animal experiments were duly approved by the Animal Ethics Committee of the National Institute of Immunology, New Delhi (IAEC/AQ/2019/185, serial no. IAEC#454/17).

For in vivo experiments, 6 to 10 weeks inbred BALB/c mice were used. They were grouped as control and infected. In the infected group, stationary phase Leishmania promastigotes were introduced through intra-cardiac injection with 50 µl 1× PBS containing 1 × 10⁶ parasites on day 0 of the study (Loeuillet et al., 2016). On day 14, the same formulation was used to give intra-peritoneal injection. At the end of day 21, the infected group was sacrificed. Another group was made from the infected BALB/c mice that were treated with miltefosine after the infection was maintained, i.e. on day 21 of infection. Five daily doses of 3 mg/kg was given orally (Dorlo et al., 2012). At the end of day 26, the treated mice were sacrificed.

Detection of L. donovani DNA in mouse tissue was carried as described elsewhere (Nicolas et al., 2002). For this, DNA was isolated from the tissues (25 mg) using DNeasy tissue kit. The following forward and reverse primers, 5’CCTATTTTACACCAACCCCCAGT-3 and 5’-GGGTAGGGGCGTTC TGCGAAA-3’ respectively, were used to amplify a 120-bp fragment of the minicircle kDNA of L. donovani; 10,000 copies of which are present in each parasite. These primers match the conserved sequences of the kinetoplast minicircle but do not match the mouse frequent nucleic acid sequences. A real-time hot-start PCR was performed with the LC FastStart DNA Master SYBR Green I Kit (Roche Diagnostics, Meylan, France) in an LC (Roche Diagnostics). The 12 µl reaction mixture contained 1 µl LC FastStart DNA Master SYBR Green I, 2 mM MgCl2, 10 µM each of the primers, and 100 ng of the template. Time and temperatures for PCR configuration were as follows: denaturation was done for 8 s at 95°C, amplification of 40 cycles was done at 95°C for 10 s, followed by 72°C for 8 s. The melting cycle was set at 95°C for 10 s, 67°C for 30 s and 95°C for 10 s followed by cooling at 40°C for 60 s. Slope (°C/s) was kept constant at 20 except in the last step of melting where the slope was 0.1. For fluorescence signal acquisition, channel F1 was used and the gain was set at 5. For normalization of fluorescent data, the F1/1 ratio was applied.

Gas–Chromatography was performed with GCMS-QP2010 Ultra (Make—Shimadzu) using Rxi-5 siloxane ms (dimensions: 30 m × 0.25 mm 1 d × 0.25 µm) column for separation. MS conditions were as follows: Ion source temperature was 230.0°C, interface temperature was 270.0°C, solvent cut time was 6.50 min, detector gain mode was kept relative, +0.00 kV was gained by detector whereas the threshold was 1,000. The program start time was 7 min and the end time was 50 min. The acquisition mode was set to scan mode. The event time was 0.20 s and the scan speed was 3,333. The m/z start value was 40 whereas the end value was 650.00.

All GC–MS raw data files were first converted to CDF format with the help of GC–MS post run analysis software which is complementary to the instrument for basic processing of data. These files were then individually uploaded and analyzed in XCMS Online software (www.xcmsonline.scripps.edu). Within the software, a single job was created using the GC/single quad centwave setting where the polarity was set to negative. After the job was completed, the analysed result was available for download in a zip folder format. The following data columns—mass/charge (‘mz’), retention time ‘rt’ and original intensity (‘into’) were copied from the provided excel sheet and saved in a.txt file, which was then converted to a.csv file for further processing and analysis of data in the MetaboAnalyst software (metaboanalyst.ca/faces/home.xhtml). MetaboAnalyst schematics include data normalization, followed by statistical analysis, ‘peak to pathway’ and enrichment analysis (Chong et al., 2018).

For peak picking and alignment, mass tolerance (m/z) was set to 0.25 and retention time tolerance to 5 s. Peaks of the same group were summed if they belonged to one sample. Peaks appearing in less than half of all samples in each group were ignored. The aligned peaks were reorganized into a single data matrix with the samples in rows and the variables (peaks) in columns for further analysis. Data normalization was done to reduce any systematic bias and to improve overall data consistency so that meaningful biological comparisons can be made. A combination of row-wise normalization by a pooled or average sample of the control group and auto data scaling was used for data normalization. This was observed to be the most suitable method that helped achieve a Gaussian curve, required for further analysis (Chong et al., 2018).

T-test/ANOVA was employed for 2-group/3-group statistical analyses respectively, as the data sets tend to follow a normal distribution and could possibly have unknown variances. Important features were selected by setting the threshold to 0.05. The absolute value changes between both groups are compared with fold change, therefore, the data before column normalization was used instead. Important features were selected by fold-change analysis with threshold ≥2. Statistical significance (P-value) versus magnitude of change (fold change) is represented by a volcano plot, which is essentially a variant of the scatter plot. Features having a fold change threshold of 2 (x-axis) and t-test threshold of (y-axis) 0.1 were considered for further analysis (Chong et al., 2018).

Using the peaks to pathways module in MetaboAnalyst, we were able to identify significant metabolites. The algorithm type was set to Metabolite Set Enrichment Analysis (MSEA) and the used pathway library was Mus musculus (KEGG). The significant metabolites identified were then proceeded for enrichment analysis using Over Representation Analysis (ORA), done in reference to HMDB using the Small Molecule Pathway Database (SMPDB). The fold enrichment factor of a pathway used here for comparison was calculated as the ratio between the number of significant pathway hits and the expected number of compound hits within the pathway (Chong et al., 2018).

Direction targeting Spearman rank correlation was used as a distance measure, that calculates the ranks of individual features rather than their actual values. A predefined profile of 1-2-1 was used to target features that showed anomaly during infection and remained more or less constant in control and treated samples. For this, all three groups of datasets—control, infected and treated samples were used for comparison (Chong et al., 2018).

To obtain a data overview, multivariate analysis was done to acquire a principal component analysis (PCA) plot, where the correlation among the datasets was obtained through a two-dimensional (2D) graph ( Figures 3A–G ). The 2D score plot draws a 95% confidence region for each group. The principal components were scaled to view separation among the control, infected and treated datasets of different tissue and biofluid samples. Results showed significant variations among the groups, which was clearly visible even in different tissue types and biofluids. A stark difference between the peak intensity patterns of the control, infected and treated groups was evident from the heat map ( Figures S1A–G ), further corroborating the PLS-DA ( Figure 3 ) and suggesting that these are indeed three separate groups with significant variations.

In case of brain samples, a fold change analysis indicated that infected samples were downregulated in majority of the features, and the degree of downregulation was much higher than that of the upregulated features in comparison to control samples. A total of 145 significant features were identified for the brain samples via volcano plot, where 88 features were upregulated and 91 features were downregulated ( Figure S2A ), whereas, in treated brain samples two upregulated and seven downregulated features were noticed as compared to control samples ( Figure S3A ). The downregulation of the features in infected brain samples was also noticeable in the heatmap ( Figure S1A ). More than 100 significant features were identified in liver samples. A clear demarcation was also visible in liver samples ( Figure 3B ), where a total of 102 upregulated and 10 downregulated significant features were identified while comparing infected vs control samples ( Figure S2B ). On the other hand, treated liver samples showed only four upregulated and 11 downregulated features in the volcano plot ( Figure S3B ). Infected splenic samples showed higher intensity features as compared to the control ones; the difference is discernable in the heat map ( Figure S1C ). A little more than 100 significant features were identified that had a threshold greater than 0.05. Majority of the infected group features seem to be upregulated with respect to the control groups in splenic samples, where 62 upregulated features and 47 downregulated features were identified ( Figure S2C ). A similar pattern was also visible in treated splenic samples where 13 downregulated features and a single upregulated feature were identified against control ( Figure S3C ). A clear separation of all the groups was observed with 2D ( Figure 3D , Figure S1D ) and 3D pls-da (figure not shown) in bone marrow samples. Three downregulated and seven upregulated features were identified on comparison of control and infected samples ( Figure S2D ). In the case of control vs treated bone marrow samples, two upregulated and five downregulated features were observed ( Figure S3D ).

Metabolic alterations were noticeable in the case of biofluids, where all of them (serum, urine and feces) showed differential and significant features. On comparison of control and infected groups, seven significant features were identified with the help of a t-test as the overall peak groups and individual features identified in serum samples were significantly less. Approximately 22 features were found to be upregulated while 18 were downregulated in the volcano plot ( Figure S2E ). Whereas, among the control vs treated group only six downregulated and four upregulated features were identified against the control ones ( Figure S3E ). In urine samples, fold change analysis suggested that the features from the infected samples were upregulated more than the control samples ( Figure S1F ). Less than 70 significant features were identified with the statistical t-test analysis. The volcano plot showed 83 upregulated and 44 downregulated features as compared to control ( Figure S2F ). Whereas in the treated vs control urine samples only a single downregulated feature was identified ( Figure S3F ). In fecal samples, two completely different sets of patterns seem to emerge from the control and infected groups ( Figure S1G ). Approximately 72 features were identified to be significant with t-test. The volcano plot shows that approximately 29 features were upregulated whereas 68 features were downregulated suggesting that majority of the significant features selected were downregulated ( Figure S2G ). The fecal samples of the treated groups showed all six significant features in the downregulated category ( Figures S3G ).

Metabolite set enrichment was done to identify and interpret patterns of metabolite concentration changes, which will further support the identified metabolites in absence of MS–MS quantification. Results showed maximum alterations in brain, urine and fecal samples during infection with respect to control ( Figure 4 ). The bone marrow of infected samples showed almost no significant changes in the metabolite sets indicating a minimum alteration in their features during infection. In the case of infected versus treated serum and spleen samples, no significant changes were noticed in the metabolite sets which indicates VL-induced changes persist for a long time even after treatment in both sources. On the other hand, in control versus treated serum and urine samples, no significant enriched metabolite-sets were noticed showing that these biofluids quickly respond to treatment and normalize back to the control configuration. Most of the pathways showed less to no enrichment in the control vs treated category showing the progressive normalization process. The pathways that remain enriched in all tissue samples even after treatment are D-glutamine and D-glutamate metabolism, arginine biosynthesis, glyoxylate–dicarboxylate metabolism, valine–leucine–isoleucine biosynthesis and histidine metabolism. This indicates that a long time is taken by the host system to normalize these pathways even after minimizing the parasitic load.

Significantly identified metabolites of MSEA from the control vs infected group were considered for pathway enrichment by applying over representative analysis (ORA) as focussing on pathway enrichment during infection is crucial to identify the systemic changes that occur during VL progression.

The identified results were further filtered through HMDB (Human Metabolome database) to get a more vivid picture of modified pathways during VL. Metabolite abundance was found to be altered for multiple central metabolic pathways including amino acids, glycolysis, TCA, fatty acid synthesis and several other pathways ( Figure 5 ). More than a few amino acid pathways like glycine, serine, methionine, arginine, proline, phenylalanine, tyrosine were found to be enriched in all tissue and biofluid samples. D-Arginine and D-Ornithine metabolism were highly altered in approximately all the tissue and biofluid samples except in serum, where no changes were perceived at all. Liver and spleen showed the maximum number of amino acid altered pathways ( Figure 5A ). Among energy-related pathways, the malate-aspartate shuttle was found to be most enriched in liver and spleen samples ( Figure 5B ). Glycolysis was also affected in all sample sources in a similar fashion, except in serum, where it was highly enriched ( Figure 5B ). De-novo triacylglycerol biosynthesis, cardiolipin biosynthesis and glycerolipid metabolic pathways followed a similar pattern for all samples. In serum and brain, no significant changes were detected in fatty-acid related metabolic pathways ( Figure 5C ). Urea cycle was found to be a highly altered pathway in liver, spleen and faeces ( Figure 5D ). Urine samples showed lesser changes in folate, porphyrin, glutathione, nicotinamide and urea cycle as compared to control.

Tissue-specific changes during Leishmania infection are an important aspect of comprehending the VL’s systemic effects. For this, several tissue-specific unique metabolites ( Tables 1 – 3 ) and common metabolites ( Figure 6 ) across systems were identified that were significantly altered in the infected samples as compared to control ones. Some of these metabolites remained altered even after treatment, showing the tendency of the host towards the disease environment. Therefore, this category of metabolites is very important in understanding the disease after-effects that need to be considered to address the VL systemic repercussions. The majority of metabolites decreased during infection as shown in Figure 6 . Only a few metabolites like glutamic acid and O-acetylserine showed an elevated pattern in serum and feces of infected BALB/c mice. Other molecules that were in higher amounts in infected mice’s brain samples are hydroxypyruvic acid, malonic acid, tartronate semialdehyde, acetoacetic acid, succinic acid semialdehyde and methyl-3-oxopropanoic acid. Liver and spleen mostly showed declining features after infection ( Figure 6 ). The highest number of unique metabolite alterations were observed in the brain ( Table 1 ), followed by spleen ( Table 2 ) and urine ( Table 3 ). Very few metabolites were identified in liver ( Table 1 ), serum ( Table 2 ) and feces ( Table 3 ) that were specifically unique to these tissues. No unique modifications were noticed in bone marrow samples. A higher rate of downregulation in metabolites of the infected samples point towards the suppressive behaviour of VL.

To understand the disease progression and modulation in host, we looked for the tissue-specific changes while witnessing the control, infected and treated samples altogether. For this, pattern hunter analysis was done using the spearman rank correlation algorithm, which was set to find a 1-2-1 pattern. This led to the identification of several metabolites, that were specifically modified in infected samples.

After an analysis of brain samples with pattern hunter, arginine, fumaric acid, oxalosuccininc acid and phosphopyruvic acid were found to be upregulated specifically during infection, whereas 4,5-Dihydroorotic acid and selenocystathionine significantly declined during infection. These metabolites were found to regain a pattern similar to that of control samples after treatment ( Figure 7A ). In liver samples, choline dramatically increased during infection whereas creatinine and oxalosuccininc acid decreased during the duration of infection. These metabolites had comparable normalized intensity in control and treated samples ( Figure 7B ). In splenic samples, 9(10)-EpoME was the only metabolite that showed increment during infection while others like glutamic acid, creatinine and dIMP decreased during the infection ( Figure 7C ). During the course of infection, N-acetylputrescine and 9(10)-EpOME were the only ones found to be upregulated in bone marrow samples that returned to control configuration after treatment ( Figure 7D ).

The major objective behind profiling biofluids was to find leads for possible diagnostic markers. Following this direction, several molecules were identified that were either serum, urine or feces specific. In serum samples, a marked increase in the intensity of glutamine, glutamic acid and pyridoxamine phosphate was noticed during infection whereas octanoic acid decreased during infection. Treated samples followed a pattern similar to that of control samples, with the exception of octanoic acid, which remained downregulated even after declined parasitic burden ( Figure 7E ). In urine, the maximum number of metabolites identified were modulated specifically during infection. Among them, hydrazine, deoxyribose, indoleacetaldehyde, guanidoacetic acid, oxalosuccinic acid, thymidine, nicotinate D-ribonucleotide and xanthosine showed elevation during infection in urine samples only. Although increasing, phosphorylcholine and N-acetylglutamic acid showed lower intensity during infection as compared to control and treated samples. 2,3-Diphosphoglyceric acid showed a significant decrease in infected and control samples ( Figure 7F ). In fecal samples, only two metabolites i.e. 1,3-diaminopropane and 3-methyl-2-oxopentanoic acid were uniquely identified, exhibiting increased intensity during infection ( Figure 7G ).