Comparative genomic analysis of L. donovani lesion-isolated amastigotes and derived promastigotes.

(A) Schematic overview of the parasite samples used for the various systems level analyses presented in this study. A more detailed overview is shown in Figure S1. (B - D) Genomic analyses. (B) Heatmap showing the somy score of three individual differentiation experiments using hamster-isolated amastigotes (ama1, ama2 and ama3) and their corresponding, culture-derived promastigotes analysed at passage 2 (pro1, pro2 and pro3). Samples and chromosomes are indicated on the x- and y-axis, respectively. The somy scores were calculated as described in Material and Methods and correlate to the grey level according to the shown legend. (C) Ratio plot showing the gene coverage ratio ama versus pro (y-axis) for all genes across the 36 chromosomes (x-axis) as calculated based on median read depth normalized by the somy score. Each color represents one individual differentiation experiment using ama isolates from individual hamsters. The red dotted lines indicate the ratio values corresponding to 1.5 (upper line) and 0.5 (lower line), indicating respectively gain or loss of one gene copy in the ama samples. Fluctuations between 0.5 and 1.5 were not considered significant. (D) Correlation plots representing the individual (histograms on the diagonal) and pairwise (off-diagonal scatterplots) distributions of SNP frequencies for the indicated samples. Of note, only SNPs with a frequency above 10% were plotted. Due to bottle neck events, some low frequency SNPs appear to be unique in one or the other stage as they are either ‘lost’ (filtered out) or ‘gained’ (passing the 10% cutoff) between the ama and pro samples. The absence of SNPs at 100% is explained by the heterozygous nature of the Ld1S genome.

Stage-specific transcript profiling.

(A) Volcano plots (left panels) showing differential transcripts abundance between ama and pro samples as assessed by RNA-seq analysis. The dotted lines indicate fold change (FC) = 2 (vertical line) and p-value = 0.01 (horizontal line). Transcripts with FC < 2 or adjusted p-value > 0.01 are represented by black dots. Transcripts with significant increased abundance FC ≥ 2 and adjusted p-value < 0.01 in ama and pro are indicated respectively in dark cyan and dark grey for all transcripts (left panel) or transcripts of only coding genes (middle panel). Transcripts for non-coding genes with significant increased abundance FC ≥ 2 and adjusted p-value < 0.01 in ama and pro are indicated in the right volcano plot (orange dots). The pie chart (right panel) depicts the percentage of transcripts for each of the categories. (B) GO term enrichment analysis of these transcripts for the category ‘biological process’ was performed, and ‘cluster efficiency’ (blue) and ‘enrichment score’ (orange) were plotted for pro (left panel) and ama datasets (all transcripts, middle panel; transcripts for coding genes, right panel). The graphs show some of the most significantly enriched GO terms associated with a Benjamini and Hochberg (BH) p-value < 0.05 that were identified using the BiNGO plug-in in Cytoscape. (C) Volcano plot corresponding to the differential expression of snoRNA between pro and ama stages. Cyan dots (ama) and grey dots (pro) represent signals with FC ≥ 2 and adj. p-value < 0.05. Black dots represent non-significant changes. (D) Localization of differentially modified rRNA sites on the Leishmania ribosome. The RNA modifications are indicated as space filling in the Leishmania cryo-EM structure (PDB: 8RXH). The identity of the RNA modification is defined by the label (Ψ, pseudouridine; Gm and Cm refer to guanine and cytosine methylation, respectively; numbers indicate residues). Upregulated and downregulated RNA modifications in amastigotes according to Rajan et al (34) are indicated in red and blue, respectively.

Systems analysis of the stage-specific transcriptome and proteome.

(A) Volcano plot showing differential protein abundance between four ama and pro samples as assessed by label-free, quantitative proteomics analysis. The dotted lines indicate fold change (FC) = 2 (vertical line) and false discovery rate (FDR) = 0.01 (horizontal line). Proteins with FC < 2 or FDR > 0.01 are represented by light grey dots. Proteins that were reproducibly detected in only one of the two stages (considered unique) are represented by the lateral histograms plotting their relative abundance according the iBAQ value. These proteins together with proteins showing differential abundance of FC ≥ 2 and adjusted p-value < 0.01 are indicated for ama in dark cyan and for pro in dark grey (see Table 6). GO term enrichment analysis for the category ‘biological process’ was performed, and the resulting values for ‘cluster efficiency’ (blue) and ‘enrichment score’ (orange) were plotted for ama (middle panel) and pro (right panel). The graphs show the main GO terms identified using the BiNGO plug-in in Cytoscape and associated with a BH p-value < 0.05. (B) Double ratio plot showing the log2 fold change between ama and pro in transcript (x-axis) and protein (y-axis) abundances. Red dots correspond to changes in both transcript and protein abundance defined either by significant differential abundance at both levels (p-value < 0.01), or by significant RNA changes (p-value < 0.01) corresponding to proteins detected in only one of the two stages. The dotted lines indicate fold changes (FC) = 2. The infinite symbol (∞) represents proteins that are unique in the ama (upper signals) or pro (lower signals) datasets. Word cloud enrichment performed with the L. donovani LdBPK orthologs is presented for each of the 4 quadrants, Q1 to Q4 and a detailed GO enrichment analysis is given in Table 7. The font size is proportional to the number of genes per GO term and the grey scale refers to the p-value calculated for each GO term. (C) KEGG map for ribosome biogenesis. Each gene associated with expression changes indicated in Q1 to Q4 (see Figure 2C) and annotated with the GO term ‘ribosome biogenesis’ has been projected on the ribosome biogenesis KEGG map. Each gene is represented by a square with the left segment showing fold changes between ama and pro at transcript and the right segment showing fold changes at the protein level, with the differential abundance indicated by the color and its intensity (red, increase; green, decrease). Only genes quantified in both transcriptome and proteome analyses and associated with an adjusted p-value < 0.01 were considered.

Analysis of stage-specific proteasomal protein turn-over.

(A) Overview of experimental workflow. To assess early and late steps of differentiation (left part of the graph in black), comparative, quantitative proteomics analyses were performed on tissue-isolated amastigotes (ama), amastigotes at 18h into the transition to promastigotes in culture (ama-18h) and fully differentiated promastigotes after 1 passage in vitro (pro). To gain insight into stage-specific and constitutive protein degradation (right part of the graph in red), amastigotes and promastigotes were treated with the irreversible inhibitor lactacystin (ama-lacta and pro-lacta). (B) Microscopic images of amastigotes at 18h in absence of lactacystin treatment (amastigotes, left), after 48h in absence (-, middle) or presence of the inhibitor (+, right). (C) Western blot analysis. Protein extracts were obtained from 2x106 parasites after 3, 6, 18 and 48h of differentiation from amastigotes to promastigotes in presence (+) or absence (-) of lactacystin. Extracts were separated by electrophoresis, transferred onto PVDF membrane and the presence of PFR2 was revealed using an anti-PFR2 antibody (upper panel). A Coomassie blue stain of the same gel is shown as loading control (lower panel). (D) Double ratio plot comparing the log2 fold changes in protein abundance between ama-18h/ama-lacta (x-axis) and ama/ama-lacta (y-axis). Dashed lines indicate FC = 2 and the infinite symbol (∞) represent proteins that are unique to one or the other comparison. (E) Projection of the proteins stabilized in presence of the inhibitor in ama (ama-lacta, orange) and pro (pro-lacta, pink) onto the double ratio plot shown in Figure 3B. Note that ama-specific proteins are rescued from degradation in pro-lacta, while pro-specific proteins are rescued from degradation in ama-lacta. (F) List of the protein kinases stabilized after lactacystin treatment in pro or ama (columns lacta rescued) and either not detected in the total proteome (column ama vs pro, grey cell) or specifically quantified at the ama (red cell) or the pro (green cell) stage. (G) Venn diagram showing the number of proteins specifically stabilized in ama or pro or stabilized in both stages after lactacystin treatment.

Stage-specific phospho-proteomic profiling.

(A) Volcano plots corresponding to the total proteome (right panel) and phosphoproteome (left panel) analyses. Proteins and phosphosites (P-sites) only identified in one stage are presented at each side of the volcano plots. Cyan dots (ama) and grey dots (pro) represent signals with FC ≥ 2 and FDR < 1%. Light grey dots represent non-significant changes. (B) Relative phosphorylation change normalized to protein abundance. Ratio plot comparing the log2 fold changes in total protein abundance between ama and pro (x-axis) and phosphosite abundance in ama versus pro (y-axis). Dashed lines indicate a FC = 2. The color intensity of each dot reflects the p-value calculated for the relative phosphorylation change normalized to protein abundance as indicated in the graph. Confidence values were derived as described in Supplementary Methods. (C) GO term enrichment analysis for the category ‘biological process. Only stage-specific phosphosites were considered (i.e. sites that showed a significant increase in relative phosphorylation normalized to protein abundance, and sites that were only detected in one or the other stage, whether the protein was identified or not in the total proteome analysis). The histogram plots show ‘cluster efficiency’ in blue and ‘enrichment score’ in orange for the GO term enrichment analysis performed with the ama (middle panel) and pro (right panel) datasets. The main GO terms identified using the BiNGO plug-in in Cytoscape and associated with a BH p-value < 0.05 are shown. (D) Phosphorylation pattern of MAPK and MAPKK proteins (left panel) or Ubiquitin hydrolases and transferases (right panel) identified in the relative phosphorylation change normalized to protein abundance. The positions of the phosphorylated amino acid specific to ama (red) and pro (green) are indicated. Only phosphopeptides identified in the relative phosphorylation change normalized to phosphorylation level (see Table 11) and detected either in one stage or the other or quantified with a FC ≥ 2 were considered for the identification of the phosphosites.

Network analyses.

(A) Model of Leishmania gene expression regulation during stage differentiation. (B-C) Network analysis. Networks, restricted to the protein kinases (blue nodes), proteins implicated in ribosome biogenesis and ribosomal proteins (violet nodes) and proteins linked to ubiquitination/deubiquitination (orange nodes) identified based on normalized phosphorylation levels conducted in ama (B) or pro (C), were generated with the STRING plug-in of the Cytoscape software package using L. infantum orthologs and full STRING network with a confidence score cutoff of 0.4. Each node represents a phosphoprotein with the respective gene identifier indicated. GO terms associated to the proteins are represented by the colored segments around the nodes according to the legend shown in the graph. Only GO terms associated with a p-value < 0.05 were considered.

Overview of the samples used in this study.

Each L. donovani infected hamster is identified by the cage number, as are the hamster-derived amastigotes (ama) and corresponding promastigotes (pro).

(A) Box plots representing the median somy score for all the chromosomes. Each box corresponds to one of three biological replicates for ama and pro as indicated in the graph. Sample identifiers are detailed in the legend of Figure 1A and Figure S1. (B) Frequency distribution plots representing, for each biological replicate and each chromosome, a kernel density estimate of the distribution of SNP frequencies in the amastigote and promastigote sample.

(A) Hierarchical clustering of the samples. The dendogram is built using the Ward method. The normed data was log-transformed first and at max 5,000 most variable features were selected. (B) Principal component analysis. The two main components are represented. Ama and pro samples correspond to blue and orange dots respectively.

(A) Cluster analysis of metabolites showing differential abundance in ama and pro samples. (B) Quantification of amino-acids or metabolites involved in glycolysis, TCA cycle and mannogen cycle. The y-axis represents the intensity for each metabolite as detected by MS. (C-D) Correlation between metabolomic and proteomic analyses for TCA cycle (C) and Glycolysis/Gluconeogenesis (D). Fold changes in proteins abundance (squares) and metabolite abundance (circles) with adjusted p-values of respectively < 0.01 and < 0.05 were projected on the respective KEGG maps. The color intensity reflects the fold change value as indicated by the legends in the graphs. (E) Representation of the mannogen pathway. Metabolites are represented by a square, with red-lined squares corresponding to the metabolites that accumulate at the ama stage. Proteins are represented by a circle, with fold changes in protein abundance between ama and pro datasets indicated by color and intensity (red, increase; green, decrease).

Label free quantitative proteomic analyses in presence or absence of lactacystin.

(A) to (E) Volcano plots representing changes in protein abundance for the different comparisons (left panels) and GO term enrichment analyses for the indicated sample (right panels). Proteins identified by at least two peptides in at least three out of four biological replicates were considered. For volcano plots, colored dots indicate values with FDR < 0.01 and FC ≥ 2 (see Table 8). The light grey dots indicate non-significant expression changes. The colored bars indicate unique protein identifications in one (green) or the other (orange) condition with relative abundance indicated by the iBAQ value (left panels). Proteins quantified in only one condition or showing a significant increase in abundance (FC ≥ 2 and adjusted p-value < 0.01) were used to perform the GO analysis for the category ‘biological process’. The histogram plots show ‘cluster efficiency’ in blue and ‘enrichment score’ in orange (right panels) for a selection of GO terms identified using the BiNGO plug-in in Cytoscape and associated with a BH p-value < 0.05. (A) ama vs ama-18h; (B) ama-18h vs pro; (C) ama vs ama lacta; (D) ama-18h vs ama lacta; (E) pro vs pro lacta.

(A) Functional network enrichment analysis for proteins showing increased abundance in ama-18h compared to ama. The network was generated with the L. infantum orthologs using the STRING plug-in of the Cytoscape software package and considering the full STRING network with a confidence score cutoff of 0.4. The GO terms associated with each protein are represented by the colored segments according to the legend. Only GO terms associated with a p-value < 0.05 were considered. (B) KEGG map for ribosome biogenesis showing changes in abundance protein between ama and ama-18h. Only proteins with an adjusted p-value < 0.01 were considered. Differential abundance is indicated by the color and its intensity (red, increase; green, decrease).

(A) Viability assessment of the indicated amastigote (upper panel) and promastigote (lower panel) samples after lactacystin treatment by FACS analysis using propidium iodide or YoPro staining. (B) Microscopic images of promastigotes at 18h in absence (upper image) or presence of 10 µM of lactacystin (lower image). (C) Western blot analysis of protein extracts obtained from ama (upper image) and pro (lower image) in presence or absence of lactacystin treatment. 10 µg of proteins were labelled with Cy5, separated on NuPAGE 4-12%, and transferred onto PVDF membrane. Membranes were subsequently incubated with a monoclonal antibody against ubiquitin and revealed with HRP-conjugated secondary antibodies. Images of Cy5-labelled proteins (left images), ubiquitinated proteins (middle images) and an overlay of both images (right images, Cy5-labelled proteins in green, ubiquitinated proteins in red) are presented. The table shows the relative quantification after normalization using the Cy5-labelled proteins of untreated samples as a loading control. (D) Western blot analysis of protein extracts obtained at the indicated time points after in vitro differentiation from hamster-derived amastigotes to promastigotes. 10 µg of protein were separated on NuPAGE 4-12%, transferred onto PVDF membrane subsequently incubated with antibodies against paraflagellar rod proteins 1 and 2 and alpha-tubulin, and revealed with HRP-conjugated secondary antibodies (left image). Micrographs showing hamster-derived amastigotes during in vitro differentiation to promastigotes (right image). Parasites were collected at the time points indicated, fixed in 4% PFA, and seeded on poly-L-lysine-treated coverslips. Parasites were observed with a 63x oil immersion objective and images were processed with AxioVision Rel.4.8 software. Differential interference contrast (DIC) microscopy was performed with Axioplan 2 imaging microscope using a 63x oil immersion objective, Axiovision software and AxioCam MRm camera (Carl Zeiss). The scale bars correspond to 1 µm.

Functional network enrichment analysis for the proteins stabilized in presence of lactacystin in promastigote (pro) and amastigote (ama) datasets.

The network was generated with the L. infantum orthologs using the STRING plug-in of the Cytoscape software package and considering the full STRING network with a confidence score cutoff of 0.4. The GO terms associated with each protein are represented by the coloured segments according to the legend. Only GO terms associated with a p-value < 0.05 were considered.

List of protein kinases stabilized in presence of lactacystin.

Their corresponding orthologs in L. donovani strain LdBPK and L. mexicana are indicated, as well as the knock-out phenotype as published by Baker et al (2021) (102).

Upset plot representing the number of phosphopeptides (y-axis, blue histograms, numbers are indicated) and the samples in which they were quantified (x-axis).

Pie chart of the phosphopeptide analysis showing the distribution of the phosphorylation sites according to the condition in which they were quantified as indicated in the legend of the graph.

Functional Network analysis.

Networks for proteins that show unique or increased phosphorylation in ama (upper panel) and in pro (lower panel) datasets were generated with the STRING plug-in of the Cytoscape software package using L. infantum orthologs and considering the full STRING network with a confidence score cutoff of 0.4. The nodes represent phosphoproteins associated with one or more GO terms (listed in the graph), with the respective gene identifier indicated. The GO terms associated to the proteins are represented by the colored segments around the nodes according to the legend shown in the graph. Only GO terms associated with a FDR < 0.05 were considered.