Introduction

During their infectious cycle, many protist and fungal pathogens undergo complex developmental transitions that adapt their biology to different host environments, thus ensuring transmission, persistence or immune evasion. For example, stage differentiation of the malarial parasite Plasmodium spp is tightly controlled at transcriptional levels by the dynamic interplay between epigenetic regulators that remodel the chromatin landscape to accommodate stage-specific transcription factors, such as AP2 DNA-binding family members (1, 2). Likewise, the transition from fast-growing tachyzoites to slow-growing bradyzoites during mammalian Toxoplasma infection is regulated by the single master TF BFD1 (3). Unlike these apicomplexan parasites, stage-differentiation in the trypanosomatid pathogen Leishmania spp is not regulated at the level of transcriptional control, raising the question of how these parasites establish and maintain the various adaptive forms that develop inside their insect and mammalian hosts.

Leishmania spp are protist parasites of humans that represent a global public health problem causing a series of immunopathologies termed the leishmaniases (4–7). These parasites differentiate into various developmental forms that are adapted for extracellular proliferation inside the sand fly midgut (procyclic promastigotes), transmission from the vector to the mammalian host during uptake of a blood meal (metacyclic promastigotes), and intracellular proliferation inside fully acidified macrophage phagolysosomes (amastigotes) (8). Unlike other eukaryotes, Leishmania development is not controlled transcriptionally as protein-coding genes in these parasites lack individual promoters and are constitutively transcribed at all stages involving parasite-specific processes such as polycistronic transcription and trans-splicing (9, 10). Despite the largely constitutive gene transcription, different Leishmania stages are characterized by striking morphological, biochemical and metabolic changes (11–16). These are regulated by a plethora of post-transcriptional mechanisms, including differential mRNA turn over through the binding of various proteins and protein complexes to the 5’ and 3’ untranslated regions (UTRs) (17–22), that can trigger mRNA degradation via deadenylation and decapping (23, 24). Differential gene expression in Leishmania can further be regulated at the translational level, for example through (i) selective ribosome recruitment by of one of the six isoforms of the Leishmania cap-binding protein eIF4E that bind to the 5’ UTR, (ii) stabilizing the poly(A) tail via the recruitment of Poly(A) Binding Protein 1 (PABP1) to the 3’ UTR promoting translation, or (iii) the differential expression of ribosomal proteins that can increase the landscape of specialized ribosomes (25–34). Finally, differential gene expression in Leishmania can also be controlled at the level of protein stability, through stage-specific expression of proteasomal components, autophagy-related genes (ATGs) or lysosomal proteases (35–41).

While each level of gene expression control has been studied individually in different Leishmania species, how these different types of regulation are integrated during Leishmania stage differentiation remains to be elucidated. Here, we approached this open question by applying a 5-layer, systems level analysis on the two major life cycle stages of Leishmania donovani, i.e. hamster-isolated, bona-fide amastigotes and derived promastigotes. Our analysis draws a complex picture of the Leishmania differentiation process that emerges from complex, co-regulated genetic networks involving mRNA turnover, protein translation, protein phosphorylation and protein degradation. Our data reveal Leishmania as an interesting model system to analyze phenotypic adaptation in the absence of transcriptional control and to dissociate the role of complex regulatory networks in the development of stable, disease-causing life cycle stages in eukaryotic pathogens.

Material & methods

Ethics statement

Work on animals was performed in compliance with French and European regulations on care and protection of laboratory animals (EC Directive 2010/63, French Law 2013-118, February 6th, 2013). All animal experiments were approved by the Ethics Committee and the Animal welfare body of Institut Pasteur and by the Ministère de l’Enseignement Supérieur, de la Recherche et de l’Innovation (projects n°#19683 and #240013).

Animals

Twenty-five female Golden Syrian hamsters (Mesocricetus auratus RjHan:AURA, weighting between 50 – 60 g) were purchased from Janvier Laboratories. All animals were handled under specific, pathogen-free conditions in biohazard level 3 animal facilities (A3) accredited by the French Ministry of Agriculture for performing experiments on live rodents (agreement A75-15-01).

Parasites and culture

Leishmania donovani strain 1S2D (MHOM/SD/62/1S-CL2D) was obtained from Henry Murray, Weill Cornell Medical College, New York, USA and maintained by serial passages in hamsters. Anesthetized hamsters were inoculated by intra-cardiac injection of 5x10⁷ amastigotes purified from hamster infected spleens as described previously (42). The weight of the animals was recorded over time, and the animals were euthanized by CO₂ asphyxiation before they reached the end-point of infection represented by a 20% loss of body weight. Amastigotes were then recovered from the infected hamster spleens and used for nucleic acid and protein extractions, or differentiated into promastigotes at 26°C in M199 complete medium (M199, 10% FBS, 25 mM HEPES; 100 µM adenine, 2 mM L-glutamine, 10 µg/ml folic acid, 13.7 µM hemin, 4.2 mM NaHCO₃, 1xRPMI1640 vitamins, 8 µM 6-biopterin, 100 units penicillin and 100 µg/ml streptomycin, pH 7.4). Promastigotes, derived from splenic amastigotes, were collected after 2 in vitro passages for nucleic acid and protein extractions.

Experimental design

Strains issued from independent experimental evolution assays are identified by number, e.g. amastigote ama1 and promastigote pro1 are the parasites prepared from hamster H1 (see Figures 1A and S1 for details). For comparative analyses, DNA, RNA, proteins and metabolites were extracted from different splenic ama and their matching pro parasites, the latter being collected for extraction during exponential growth phase.

Figure 1:

Parasite differentiation and lactacystin treatment

Spleen-derived amastigotes were incubated at 2x10⁷/ml and 26°C in presence or absence of 10 µM of lactacystin (L6785, Sigma) for 3, 6, 18 and 48h in M199 complete medium. At each time point, parasites were collected to control for the presence of paraflagellar rod protein 2 (PFR2) as promastigote differentiation marker. After one in vitro passage, amastigote-derived promastigotes were incubated at a concentration of 2x10⁷ parasites per ml in presence or absence of 10 µM of lactacystin for 18h. Parasite viability was assessed by FACS analysis after 18h of treatment by propidium iodide (1µg/ml) or YO-PRO-1 (0.2 μM) (# Y3603, Invitrogen) staining as readout using untreated and PFA-treated parasites as controls. Parasites were collected from 25 ml of culture after 18h of treatment for quantitative, label-free proteomics analyses.

Microscopy

Differential interference contrast (DIC) microscopy was applied on hamster-derived amastigotes during in vitro differentiation to promastigotes using an Axioplan 2 imaging microscope using a 63x oil immersion objective, the AxioVision Rel.4.8 software and an AxioCam MRm camera (Carl Zeiss). Alternatively, microscopic images of amastigotes in presence or absence of the lactacystin inhibitor were obtained using the EVOS FL microscope (life technologies) at a 20x magnification.

Western blot analysis

Protein extracts were obtained from parasites after 3, 6, 18 and 48h of differentiation from amastigotes to promastigotes in presence or absence of lactacystin. Extracts were loaded on 4-12% SDS-PAGE (NuPage 4-12% Bis-Tris gel NP0321BOX, Invitrogen), separated by electrophoresis and stained with Coomassie blue or transferred onto PVDF membrane. The presence of PFR2 and α-tubulin was revealed using the anti-PFR2 antibody (kindly given by Philippe Bastin), an anti- α-tubulin antibody (clone B-5-1-2, Sigma), the corresponding HRP-conjugated secondary antibodies (Invitrogen) and the SuperSignal West Pico PLUS kit (Thermoscientific). Protein extracts from lactacystin-treated or untreated amastigotes and promastigotes were labelled with Cy5, loaded on 4-12% SDS-PAGE (NuPage 4-12% Bis-Tris gel NP0321BOX, Invitrogen), separated by electrophoresis and transferred onto PVDF membrane. The presence of ubiquitinated proteins was revealed using the mouse anti-ubiquitin monoclonal antibody P4D1 (Abcam), the HRP-conjugated anti-mouse secondary antibody (Invitrogen) and the ECL Prime western blotting detection reagent (Cytiva). Images were acquired using the ImageQuant 800 (Cytiva).

Genome sequencing and data analysis

DNA was prepared from splenic amastigotes and promastigotes at exponential culture phase (three biological replicates). Parasites were centrifuged at 1,600 x g (pro) or 2,000 x g (ama) for 10 min at room temperature. Approximately 1x10⁸ promastigotes from logarithmic growth phase were re-suspended in 200 µl PBS and genomic DNA was purified using DNeasy Blood and Tissue kit from Qiagen and RNase A, according to the manufacturer’s instructions. DNA concentrations were measured in duplicate by fluorescence using a Molecular Device fluorescence plate reader (Quant-IT kits, Thermo Fischer Scientific). The quality of DNA was controlled by determining the DNA Integrity Number (DIN) analyzing 20 ng of DNA on a TapeStation 4200 (Agilent). One µg of genomic DNA was used to prepare a library for whole genome sequencing on an automated platform, using the Illumina “TruSeq DNA PCR-Free Library Preparation Kit”, according to the manufacturer’s instructions. After normalization and quality control, qualified libraries were sequenced on a HiSeqX5 platform from Illumina (Illumina Inc., CA, USA) at the Centre National de Recherche en Génétique Humaine (CEA, Evry, France), generating paired-ended, 150-bp reads. Sequence quality parameters were assessed throughout the sequencing run. Standard bioinformatics analysis of sequencing data was based on the Illumina pipeline to generate a FASTQ file for each sample.

Genomic DNA reads were aligned to the L. donovani Ld1S reference genome (https://www.ncbi.nlm.nih.gov/bioproject/PRJNA396645, GCA_002243465.1) with BWA mem (version 0.7.12) with the flag -M to mark shorter split hits as secondary. Samtools fixmate, sort, and index (version 1.3) were used to process the alignment files and turn them into bam format (43). RealignerTargetCreator and IndelRealigner from the GATK suite were run to homogenize indels (44). Eventually, PCR and optical duplicates were labeled with Picard MarkDuplicates (version 1.94(1484)) (https://broadinstitute.github.io/picard/) using the option “VALIDATION_STRINGENCY=LENIENT”. For each read alignment file, Samtools view (version 1.3) and BEDTools genomecov (version 2.25.0) were used to measure the sequencing depth of each nucleotide (45). Samtools was run with options “-q 50 -F 1028” to discard reads with a low map quality score or potential duplicates, while BEDTools genomecov was run with options “-d -split” to compute the coverage of each nucleotide. The coverage of each nucleotide was divided by the median genomic coverage. This normalization is done to account for library size differences. The chromosome sequencing coverage was used to evaluate aneuploidy between EP.1 and LP.1 samples. Then for each sample and for each chromosome, the median sequencing coverage was computed for contiguous windows of 2,500 bases. As previously published (46), the stably disomic chromosome 36 was used to normalize chromosome read depth and to estimate chromosome polysomy levels in each sample. Gene counts were produced using featureCounts (version 1.4.6-p3 (47)) with these parameters: -s 0 -t gene - g gene_id and were normalized according to the median-ratio method.

For the generation of the chromosome median somy score heatmap, mean coverage in 300 bp bins as generated by the GIP pipeline (48) were used to compute somy scores per chromosome by first normalizing bin scores for a sample by their median across all genome (to obtain comparable values between samples), multiplying by two (to scale somy values to the default diploid state assumed for most of the genome), and then taking the median across bins belonging to a given chromosome (see also https://gip.readthedocs.io/en/latest/giptools/karyotype.html). These chromosome scores were then displayed on a sample x chromosome heatmap, excluding the “maxi” chromosome. This was done using the Pandas 1.4.2 (49, 50) Pandas-dev/pandas: Pandas 1.2.4, https://zenodo.org/records/4681666 (50)), Matplotlib 3.5.1 (51) and Seaborn 0.11.2 (52) Python libraries.

For the assessment of gene coverage ratios, normalized mean coverages per gene are reported by GIP (version 1.1.0 (48)), as the mean coverage of the gene divided by the median coverage of the chromosome containing the gene. For a given gene, and a given amastigote - promastigote pair, a gene coverage ratio was computed adding 0.1 to the gene coverages of the individual samples, in order to avoid zero division errors. The genes were sorted by genomic coordinate and assigned successive integers to create a “genomic index”. Gene coverage ratios were then plotted along this genomic index, using the same Python libraries as above.

For the SNP analysis, alternative (Alt) allele frequencies for SNPs were retrieved from the output of the SNV giptools module (provided along with the GIP pipeline). Euclidean distances between samples in the SNP allele frequencies space were computed and used to build a dendrogram using the scipy.cluster.hierarchy Python module (SciPy version: 1.8.0). The clustering was obtained using the UPGMA method.

RNA-seq analysis

Total RNA was extracted from splenic amastigotes and promastigotes in exponential culture phase (4 biological replicates). Amastigotes and promastigotes were centrifuged respectively at 2,000 x g or 1,600 x g for 10 min at 4°C, and re-suspended in the lysis buffer supplied with the Qiagen RNeasy Plus kit. The samples were stored at -80°C and RNA extractions were performed according to the manufacturer’s instructions, including a DNase treatment. RNA integrity was validated using the Agilent Bioanalyzer. DNase-treated RNA extracts were used for library preparation. Libraries were constructed using an Illumina Stranded mRNA Prep (Illumina, USA) following the supplier’s recommendations. RNA sequencing was performed at the Biomics Center (Institut Pasteur, Paris, France) on the Illumina NextSeq 2000 platform for a target of 40M paired-end reads per sample.

The RNA-seq analysis was performed using the Sequana RNA-seq pipeline (version 0.20.1, https://github.com/sequana/sequana_rnaseq) from the Seauqnq project (53). The pipeline workflow manager is Snakemake 7.3.2 (54). The pipeline was executed with default parameters and bioinformatics software are available as reproducible containers from the Damona project (https://github.com/cokelaer/damona). Reads were mapped to the Leishmania donovani reference using bowtie2 2.4.5 (55). Gene-level quantification was conducted using FeatureCounts 2.0.1 (47) assigning reads to genomic features based on the gene annotation and gene_id attribute, while accounting for strand-specificity information. Differential expression analysis was performed using DESeq2 v1.24.0 (56) scripts available in the Sequana library. Statistical testing identified differentially expressed genes by comparing amastigote and promastigote sample groups, with significance determined using Benjamini-Hochberg adjusted p-values (false discovery rate FDR < 0.05).

To specifically assess the differential expression of snoRNA, RNA-seq libraries from promastigote and amastigote stages were mapped to the Ld1S genome, and snoRNA read counts were obtained using bedtools multiBamCov. For snoRNAs with multiple copies, read counts were summed per sample to generate a single expression value per snoRNA family. DESeq2 was then used to perform differential expression analysis between amastigote and promastigote stages based on the summed snoRNA counts.

Label-free, quantitative total proteome analyses

Untreated amastigotes (ama and ama-18h) and promastigotes (pro), and lactacystin-treated parasites (ama lacta and pro lacta) were recovered from four replicates (Figure 1A and S1) and washed three times with cold PBS at 2,000 x g or 1,600 x g for 10 min at 4°C. Parasite lysates were prepared in eFASP lysis buffer (4% SDS / 0.2 % DCA / 50mM TCEP / 50 mM ammonium bicarbonate buffer pH 8) (61). For detailed protocol see supplementary data file. MS scans were acquired at a resolution of 70,000 and MS/MS scans (fixed first mass 100 m/z) at a resolution of 17,500. The AGC target and maximum injection time for the survey scans and the MS/MS scans were set to 3 x 10⁶ for 20 ms and 10⁶ for 60 ms, respectively. An automatic selection of the 10 most intense precursor ions was activated (Top 10) with a 40 s dynamic exclusion. The isolation window was set to 1.6 m/z and normalized collision energy fixed to 28 for HCD fragmentation. We used a minimum AGC target of 10⁴ corresponding to an intensity threshold of 1.7 x 10⁵. Unassigned precursor ion charge states as well as 1, 7, 8 and >8 charged states were rejected and match was disabled.

Label-free quantitative phosphoproteome analyses

Amastigote (ama) and promastigote (pro) parasites from four biological replicates were recovered and washed three times by centrifugation in cold M199 at respectively 2,000 x g or 1,600 x g for 10 min at 4°C. Samples were incubated for 10 min at 4°C in lysis buffer (1ml per 1.5x10⁹ parasites) consisting of 8 M urea, 50 mM Tris, supplemented with a protease inhibitor cocktail (cOmplete from Roche) and a phosphatase inhibitor cocktail (PhosStop from Roche). Following sonication for 5 min using a sequence of 10 s pulse and 20 s pause, the lysates were centrifuged 15 min at 14,000 x g and 4°C, and the supernatant was collected and stored at -80°C until use. Proteins were quantified by RC DC protein assay (Bio-Rad) and adjusted to 1.3 µg.µl^-1 in lysis buffer. For detailed protocol see supplementary data file. Phosphopeptide enrichment was carried out as described in (62) and detailed in supplementary data. All analyses were performed on a Q Exactive HF Mass Spectrometer (Thermo Fisher Scientific) coupled with a Proxeon EASY-nLC 1000 (Thermo Fisher Scientific) as detailed in supplementary data. Details for the proteomics and phosphoproteomics analyses following the data acquisition are available in supplementary data. Results are presented in table 6 and 10.

Phosphopeptides were selected for Gene Ontology enrichment analyses if (i) they were exclusively identified in one of the two stages, (ii) they showed significant, stage-specific changes in phosphopeptide abundance (fold change ≥ 2, adj. p-value < 0.01) even if the corresponding protein was not detected in the total proteome analysis, and (iii) they showed a significant increase in relative phosphorylation normalized to protein abundance (p-value < 0.05) as calculated by the ratio ‘change in phosphopeptide abundance’ vs ‘change in protein abundance’ using a cut off of fold changes ≥ 2 and adj. p-value < 0.01 for both analyses (see table 11).

Metabolomic analysis

The sample extraction was performed as described previously (63). Briefly, 10⁸ cells were used per each 200 μL sample (4 replicates/stage). Firstly, cells were rapidly quenched in a dry ice/ethanol bath to 4°C, then centrifuged, washed with 1 x PBS, and resuspended in the extraction solvent (chloroform:methanol:water, 1:3:1). After vigorous shaking at 4°C for 1 h, extracts were centrifuged (16,000 x g, 4°C, 10 min) and the supernatants collected and stored at -80°C until the analysis. LC-MS analyses were performed using separation on 150 x 4.6 mm, 5 mm ZIC-pHILIC (Merck) on UltiMate 3000 RSLC (Thermo Scientific) followed by mass detection on an Orbitrap Exactive mass spectrometer (Thermo Fisher) at Glasgow Polyomics. Analyses were performed in positive and negative polarity switching mode, using 10 μl injection volume and a flow rate of 300 μl/min over 26 min on the column maintained at 30°C, as follows: 0 to 20 min 20%-to 80% solution A, 15 to 17 min 95% solution A, 17 to 26 min 20% solution A where solution A is 20 mM ammonium carbonate in water and solution B is acetonitrile). The samples were run alongside 249 authentic standards at 10 μM each. Mass spectrometry data was processed using Mzmatch (64) and Ideom (63) software. Unique signals were extracted using the centwave algorithm and matched across biological replicates based on mass-to-charge ratio and retention time. These grouped peaks were then filtered based on relative standard deviation and combined into a single file. The combined sets were then filtered on signal-to-noise score, minimum intensity and minimum detections. The final peak set was then gap-filled and converted to text for use Putative metabolite identification corresponds for the most part to Metabolite Standards Initiative (MSI) level 2 (mass only, and thus only considered a tentative annotation), whereas metabolites matching in retention time to an included standard correspond to level 1 (considered likely an accurate annotation). Peaks having an area with root squared deviation across pooled samples > 50% were excluded, as were those with a retention time < 4 min (due to poor resolution).

Localizing rRNA modification on Leishmania cryo-EM ribosome structure

The cryo-EM atomic model of Leishmania major 80S ribosomes bound to mRNA and all three tRNAs (PDB: 8RXH) was used to project rRNA modifications of the ribosome. Figures were generated using UCSF Chimera-X software (65).

Systems-level analyses

A lists of L. donovani GO terms (66) was built in house. The Biological Networks Gene Ontology tool (BiNGO) plugin of the Cytoscape software package (version 3.8.2) was used to map and visualize functional enrichments in each data set based on the GO hierarchy. A Benjamini & Hochberg false discovery rate with a significance level of 0.05 was applied. Cluster efficiency represents the percentage of genes for a given GO term compared to the total number of genes with any GO annotation in the considered set of genes. Enrichment score corresponds to the percentage of genes for a given GO term compared to all the genes sharing the same GO term in the genome. Word Cloud of GO enrichment analysis limited to GO Slim terms and a threshold of p-value < 0.05 for the category Biological Process was performed on TriTrypdb (https://tritrypdb.org) using the L. donovani LdBPK orthologs for all genes that showed statistically significant expression changes (p-value < 0.01) at both RNA or protein levels.

For the phosphoproteomic work, phosphopeptides were selected for Gene Ontology enrichment analyses if (i) they were exclusively identified in one of the two stages, (ii) they showed significant, stage-specific changes in phosphopeptide abundance (fold change ≥ 2, adj. p-value < 0.01) even if the corresponding protein was not detected in the total proteome analysis, and (iii) they showed a significant increase in relative phosphorylation normalized to protein abundance (p-value < 0.05) as calculated by the ratio ‘change in phosphopeptide abundance’ vs ‘change in protein abundance’ using a cut off of fold changes ≥ 2 and adj. p-value < 0.01 for both analyses (see table 11).

Functional enrichment networks were built with the STRING plugin of the Cytoscape software package (version 3.8.2) using the L. infantum orthologs, and the full STRING network with a confidence score cutoff of 0.4. A false discovery rate with a significance level of 0.05 was applied. The “ClusterProfiler” package (version 4.2.2) of R was used for KEGG (Kyoto encyclopedia of genes and genomes) gene set enrichment analysis (GSEA) and for data mapping on metabolic pathways available in the KEGG database. Results were visualized using the “pathview” packages (version 1.34.0).

Results

Leishmania stage differentiation occurs independently from changes in gene dosage

We carried out an in-depth, systems analysis of independent L. donovani amastigote preparations isolated from infected hamster spleens and their culture-derived promastigotes at in vitro passage 2, with the aim of revealing mechanisms controlling stage-specific differential expression across five quantifiable information levels, by comparing the genome, transcriptome, proteome, metabolome and phospho-proteome of the two life-cycle stages (Figure 1A and Figure S1). We first investigated the possible role of genome instability in stage-differentiation, given the very well documented propensity of L. donovani to establish and tolerate chromosome and gene copy number variations in culture and during sand fly or hamster infection (46, 66, 70).

Applying our computational pipeline termed ‘GIP’ (48) to the genome sequences of three independent pairs of amastigotes and promastigotes confirmed the largely disomic state of the Leishmania genome at both stages, except for the known constitutive tetrasomy of chr 31 (Figure 1B and Figure S2A, Table 1). We did not detect any major changes at karyotypic levels following promastigote differentiation, confirming our previous observations that stage differentiation occurs independent of karyotypic adaptation (46, 70). Nevertheless, a clear pattern emerged for promastigotes that showed a slight but reproducible increase in median coverage for chr 1, 2, 3, 4, 6 and 20, suggesting the emergence of mosaic aneuploidy for these chromosomes early during culture adaptation. However, these minor signals are not predictive for the karyotypic changes characteristic of culture-adapted Ld1S parasites (i.e. chr 5, 26, 33, (46, 70)). Likewise, differentiation of promastigotes was not associated with any major gene copy number variation at early passage 2 (Figure 1C, Table 2). Even though a number of genes showed increased (e.g. Ld1S_300477100 encoding an intraflagellar transport (IFT) protein) or decreased (e.g. Ld1S_220235700 encoding a hypothetical protein) read depth in individual pro strains, none of these changes converged across all pro samples. In the absence of convergent selection, it is impossible to distinguish if these gene CNVs provide some strain-specific advantage or are merely the result of random genetic drift. Plotting the normalized read depth ratios between amastigotes and promastigotes revealed an intriguing, undulating pattern that was highly reproducible across the three independent experiments and has been previously associated with nascent DNA produced during the replicative S-phase of the cell cycle (71–73). Finally, differentiation did not affect the SNP frequency distribution (Figure S2B, Table 3), which resulted in a largely diagonal pattern when plotting amastigote against promastigote SNP frequencies (Figure 1D).

In conclusion, unlike long-term culture adaptation, which drives important karyotypic changes (46, 74) that can affect transcript- and protein abundance levels as demonstrated by a recent 4-layer systems analysis (75) in vitro differentiation of splenic amastigotes into promastigotes was not associated with significant changes in chromosome or gene copy number, ruling out gene dosage effects as the source of stage-specific expression differences. Thus, L. donovani promastigote differentiation in vitro is independent of genomic adaptation, confirming our previous results obtained with sand fly-derived, promastigote parasites (66, 76).

Comparative RNA-seq analysis reveals stage-specific, co-regulated gene clusters that are controlled at post-transcriptional levels

We next analyzed stage-specific expression changes applying RNA-seq analyses on four amastigote and derived promastigote strains. In contrast to the relatively stable genome infrastructure maintained between stages, we observed substantial, stage-specific differences in transcript abundance. From 10402 detected transcripts in both stages, 6478 showed statistically significant differences (adj. p-value < 0.01, base mean reads ≥ 10 reads), including 1234 and 1027 transcripts with a 2-fold or higher increase in abundance in amastigotes and promastigotes, respectively (Table 4). Considering the lack of promoter-driven expression control of individual, protein coding genes in Leishmania and the absence of stage-specific gene dosage effects (see Figure 1), the observed differences in amastigote and promastigote transcript abundances are most likely attributable to stage-specific mRNA turnover. Interestingly, non-coding (nc) RNAs account for 64.7% of the differentially expressed transcripts at the amastigote stage (Figures 2A).

Figure 2:

Analyzing the stage-specific RNA-seq data sets identified clusters of co-regulated transcripts that share a common annotation (Figure 2B, Table 4). In promastigotes, we observed increased abundance for 66 transcripts annotated for the interconnected GO terms ‘microtubule based process’, ‘cytoskeleton organization’ and ‘ciliary cell motility’, confirming the post-transcriptional co-regulation of flagellar biosynthesis as previously shown in L. mexicana (77) (Figure 2B left panel, Table 4, sheet G). Other co-regulated, functional gene sets are defined by the GO terms ‘oxidation reduction’ (46 transcripts), ‘glycolysis’ (15 transcripts), ‘post-transcriptional regulation of gene expression’ (46 transcripts) and ‘translational elongation’ (28 transcripts), further supporting translational control as a regulatory step in stage differentiation.

In amastigotes, a major co-regulated gene cluster is represented by 77 amastin transcripts encoded on 5 different chromosomes that show increased abundance in amastigotes, confirming the stage-specific, post-transcriptional regulation of this gene family as previously described (19, 78, 79) (Table 4, sheet C, E and F). Another functional gene cluster showing increased abundance in amastigotes was associated with protein translation, with GO enrichment observed for the terms ‘ribosome biogenesis’ (26 transcripts), ‘ncRNA processing’ (22 transcripts), including 17 transcripts annotated for the GO term ‘rRNA processing’, and 6 transcripts for ‘tRNA processing’ (Figure 2B, middle and right panel, Table 4, sheet E and F).

These data suggest ribosome and epi-transcriptomic regulation as possible central components in Leishmania stage differentiation. Indeed, analyzing poly-A transcriptome libraries of amastigotes and promastigotes we observed significant differences in stage-specific expression of 110 pre-snoRNAs (Table 5, Figure 2C). These data fit with the expression of a series of methyl- and pseudouridyl-transferases at this stage (see Table 6, sheets A and E). Stage-specific expression of some of these snoRNAs correlated with stage-specific changes previously observed in Pseudouridine and 2’-O-Methylation rRNA modifications (34). The positions of these modifications and the corresponding snoRNAs known to affect ribosome translational activity are indicated on the high resolution cryo-EM structure of the Leishmania ribosome (Figure 2D).

In conclusion, our RNAseq data reveal a series of functionally related, stage-specific gene clusters whose mRNA abundance is regulated in a coordinated manner at post-transcriptional level. Whether this co-regulation is governed by stage-specific mRNA decay mechanisms acting on 5’ and 3’ UTRs or by stage-specific mRNA stabilization via RNA-binding proteins or mRNA modification remains to be established.

Quantitative proteome analysis reveals co-regulated gene clusters that are independent of stage-specific mRNA turnover

We next employed quantitative proteomic analysis to assess how the different co-regulated gene clusters observed in our RNAseq analyses translate into stage-specific proteomes. Label-free quantitative proteomic analysis of 4 replicates of amastigotes and derived promastigotes identified over 4000 proteins, including 1987 differentially expressed proteins (adjusted p-value < 0.01), and 1534 that were exclusively detected in either ama or pro (Figure 3A left panel, Table 6). Considering the 822 proteins that are specific or more abundant in amastigotes (fold change ≥ 2), functional enrichment was mainly observed for the GO terms ‘oxidation reduction’ (61 proteins), ‘vesicle mediated transport’ (23 proteins) and ‘phosphorylation’ (41 proteins) (Figure 3A, middle panel, Table 6, sheets B and D). Considering the 1637 proteins that are specific or more abundant in promastigotes (fold change ≥ 2), functional enrichment was observed for the GO terms ‘post-transcriptional regulation of gene expression’ (71), ‘translation’ (69 proteins), ‘ribosome biogenesis’ (21 proteins) and ‘ciliary cell motility’ (42 proteins) (Figure 3A, right panel, Table 6, sheet C and E).

Figure 3:

The metabolomic analysis of four independent amastigote and promastigote samples (Figure S4A) supported our proteomic data. Promastigotes showed increased amino acid biosynthetic activity that matched the changes in protein abundance observed for this pathway, while the amastigote metabolome indicates a shift from glycolytic to TCA cycle-dependent energy production that was previously described (12, 14, 80) (Figures S4B-D), with a notable increase in metabolites involved in butanoate metabolism indicative of a shift between carbohydrate use as a primary energy source to utilization of fatty acid oxidation. The amastigote expression profile and metabolome further indicate the stage-specific production of the storage-carbohydrate mannogen, thus corroborating previous observations (81) (Figure S4E).

Surprisingly, only few enriched GO terms matched between transcriptomics and proteomics datasets, suggesting additional layers of regulation between the transcriptional and translational levels of gene expression. To investigate this aspect on a genome-wide level, we next plotted the simple numerical values we obtained when calculating the ratios of expression changes between ama and pro observed in the transcriptomics and proteomics datasets (Figure 3B, Table 7). Only expression changes were considered that either showed statistically significant differential abundance at both RNA and protein levels (p < 0.01), or showed significant RNA changes (p < 0.01) with the corresponding protein being detected in only one of the two stages. These latter proteins are identified by signals that were arbitrarily placed at the upper (detected in ama) or the lower (detected in pro) parts of the graph. Whether these proteins just escape detection due to low expression or are truly not expressed remains to be established. The double ratio analysis of the resulting 2349 gene products defined four quadrants: Quadrants 2 (Q2) and 4 (Q4) correspond to 430 and 1042 genes (Table 7, sheet E and I), respectively, for which stage-specific expression changes correlate between protein and RNA levels, suggesting that the abundance of these proteins may be mainly regulated by mRNA turn-over. Q2 (mRNA and protein up in ama) includes the GO terms ‘vesicle mediated transport’ (17 genes) and ‘lipid metabolic process’ (21 genes) (Table 7, sheet F), while Q4 (mRNA and protein up in pro) includes the GO term ‘cell motility’ (41 genes) (Table 7, sheet J).

Counter-correlations between RNA and protein abundances (p-value < 0.01) were observed for 463 and 414 genes, respectively falling into quadrants Q1 (RNA down and protein up in ama) and Q3 (RNA down and protein up in pro) (Figure 3B, Table 7, sheets C and G). The most important difference was observed for 35 genes annotated for the GO term ‘ribosome biogenesis’ (Figure 3C, Table 7, sheet H). While this GO term defined a positively regulated gene cluster at mRNA level in amastigotes (see Figure 2B), the same genes show reduced expression at this stage at protein level (Figure 3C). This lack of correlation may be explained by stage-specific changes in mRNA translation efficiency.

In conclusion, while many stage-specific changes of protein abundance are strictly regulated by similar changes at mRNA level, others show discrepancies in stage-specific abundances between transcript and protein levels, revealing regulatory mechanisms that act independent of mRNA stability, for example at the level of translational control or differential protein stability.

Stage-specific protein turn-over is required for parasite differentiation

The discrepancies we observed in a sub-set of genes between changes in mRNA and protein abundance during the amastigote-to-promastigote transition (see Figure 3B) primed us to investigate the role of proteasomal protein degradation in Leishmania differentiation. We chose the highly specific and irreversible proteasome inhibitor lactacystin (82, 83) over the typanosomatid-specific, reversible drug candidate LXE408 (84) as the latter’s potent cytotoxicity can confound direct effects on protein turnover with secondary consequences of cell death, limiting its utility for dissecting proteasome function in living parasites. Comparative, quantitative proteomics analyses were performed with (i) spleen-derived amastigotes (ama), (ii) amastigotes after the first 18h of differentiation in the absence (control, ama-18h) and the presence of latacystin (ama-lacta), and (iii) promastigotes in culture without (pro) and after 18h of treatment with lactacystin (pro-lacta) (Table 8, Figure 4A and S5).

Figure 4:

This data set allowed us important insight into processes of early differentiation, revealing increased abundance in ama-18h compared to ama for 694 proteins (fold change ≥ 2, adjusted p-value < 0.01) (Table 8b, sheet B) with GO enrichment observed for the terms ‘ribosome biogenesis’, ‘RNA methylation’, ‘rRNA processing’, ‘gene expression’, and ’post-transcriptional regulation of gene expression’, the latter including 6 genes encoding for pseudouridine transferases known to modify ribosomal (r) RNA (Figures S5A and S6) (Table 8b, sheets D to H). These data further sustain a role of epitranscriptomic regulation and stage-adapted ribosomes as key processes that may initiate Leishmania promastigote differentiation, at least in culture (31, 34, 66).

Lactacystin treatment further uncovered an essential regulatory role of proteasomal protein degradation in L. donovani development, which builds on previous published observations (85). While the presence of the inhibitor neither affected parasite viability (Figure S7A) nor morphology (data not shown and Figure S7B) or promastigote growth (as previously shown by Silva-Jardim et al (85)), lactacystin treatment abrogated the amastigote-to-promastigote developmental transition as judged by the persistence of an amastigote-like morphology (i.e. oval cell shape, retracted flagellum), and the absence of PFR2 expression characteristic of promastigotes (Figures 4B and C). In contrast, untreated amastigotes spontaneously converted to flagellate promastigotes after 48h in culture (Figure S7D). Proteasome inhibition by lactacystin blocked amastigote-to-promastigote differentiation without inducing rapid global accumulation of ubiquitinated proteins (Figure S7C), which are observed only at later time points and correlated with differentiation-associated remodeling rather than acute proteasome inhibition.

Comparing the proteomics signatures of lactacystin-treated samples and untreated controls (Figure S5) allowed us to reveal the proteasomal targets, whose rescue correlated with the observed abrogation of differentiation. This analysis is complicated by the fact that amastigotes will not only respond to the 18h treatment but also initiate their differentiation into promastigotes during this period, at least in the untreated control sample. To distinguish between changes caused by lactacysin treatment and changes due to the differentiation process, we plotted the ratio of ama/ama-lacta against the ratio ama-18h/ama-lacta (Figure 4D), which defined one quadrant (Q3) representing 129 proteins whose increased abundance in ama-lacta was exclusively due to rescue from degradation (Table 9, sheet A). Applying STRING and GO analyses on these proteins, some of the most significantly enrichment terms included ‘protein kinase’, ‘cellular response to stimulus’, and ‘cell cycle’ (Figure S8 lower panel, Table 9, sheet B and G).

Proteomics analysis of lactacystin-treated promastigotes (pro-lacta) revealed the rescue of 289 proteins from degradation when compared to untreated control (pro) (Table 9, sheet C). Comparison of rescued proteins across both ama and pro stages revealed a series of intriguing findings: First, most of these proteins were rescued in a stage-specific manner, confirming differential protein degradation as an important regulatory component in the development and maintenance of Leishmania life cycle stages as suggested by previous reports (37, 86–92). Second, just like the proteins rescued from proteasomal degradation in ama-lacta, the pro-lacta data set too showed GO enrichment for the terms ‘protein kinase’ and ‘cellular response to stimulus’, involving distinct sets of proteins (Figure S8, upper panel) (Table 9, sheets D, E and G). Thus, Leishmania differentiation correlates with the expression of complex signaling networks that are established in a stage-specific manner. Third, when projecting the lactacystin-rescued proteins to the established, stage-specific proteomes (see Figure 3), an intriguing pattern emerged: the majority of proteins that were rescued in ama-lacta turned out to be part of the promastigote proteome data set (Figure 4E). Vice versa, many proteins rescued in pro-lacta were part of the amastigote-specific proteome data set (Table 9, sheets H and I), a pattern that affected 14 protein kinases (Figure 4F). Thus, these proteins seem to be expressed constitutively across both stages, while their steady-state stage-specific abundance is regulated by proteasomal degradation. Surprisingly, 12 proteins (including the protein kinase Ld1S_220237000) were rescued from degradation in both ama-lacta and pro-lacta but were never identified in the proteomes of the corresponding untreated controls, suggesting that they are constitutively degraded in both stages under our experimental conditions (Figure 4G, Table 9, sheet F).

In conclusion, our results confirm the important role of protein degradation in regulating the L. donovani amastigote and promastigote proteomes and identify protein kinases as key targets of stage-specific proteasomal activities.

Quantitative phosphoproteome analysis reveals Leishmania biological networks associated with stage differentiation

The investigations above associated the expression of protein kinases (and their degradation) to Leishmania stage differentiation (Figure 4F, S8 and S9), thus confirming previous reports on the importance of phosphotransferase activities in trypanosomatid development (93–103). Differential turnover of these important signaling proteins by stage-specific proteasomal activities likely represents a regulatory switch controlling parasite development. To further investigate how this putative regulatory switch operates, we screened for stage-specific protein kinase substrates applying label-free, quantitative phospho-proteomic analyses on three independent amastigote and promastigote samples. We identified a total of 2079 phospho-peptides in amastigotes and 7095 phospho-peptides in promastigotes, suggesting a three-fold increase in global protein phosphorylation in insect-stage parasites (Table 10). However, increased TiO₂ phosphopeptide enrichment in a given sample may be due either to increased substrate expression or increased de novo phosphorylation (Figure 5A and S10). To distinguish between these two possibilities we assessed stage-specific relative phosphorylation changes normalized to protein abundance (Figure 5B).

Figure 5:

We identified 877 phosphoproteins in amastigote samples carrying 2079 phosphosites of which 55.8% (1160) were uniquely identified at this stage, while 12.6% (565) showed a significant change in normalized phosphorylation levels compared to promastigotes (Figure S10, Table 11 sheets A to D). Conversely, 2054 phosphoproteins with 7095 phosphorylation sites were identified in promastigote samples, of which 70% (4967) were exclusively detected at this stage, while 1% (73) presented a significant change in normalized phosphorylation levels compared to amastigotes (Figure S10, Table 11 sheets F to I). These data reveal a surprising level of stage-specific phosphorylation in promastigotes, which may reflect their increased biosynthetic and proliferative activities compared to amastigotes, or may be a consequence of culture adaptation.

We next investigated the pathways regulated by phosphotransferase activities analyzing the phosphoproteomes of amastigote and promastigote for enriched biological processes (Figure 5C) (Table 11, sheets E and J). Stage-specific enrichment was observed in the pro dataset for ‘ciliary cell motility’, indicating that flagellar biogenesis and activity are not only regulated by increased protein abundance, but also by increased phosphorylation of flagellar proteins. Surprisingly, in both the ama and pro datasets, similar functions were enriched, including cell division (e.g. ‘cytokinesis’, ‘cell cycle’), signaling (e.g. ‘protein aa phosphorylation’, ‘signal transduction’), or post-transcriptional regulation of gene expression. Thus, aside stage-specific protein turnover (see above, Figure 4), differential phosphorylation represents yet another layer of regulation that contributes to the establishment of biological networks in Leishmania that share similar function in both ama and pro, but whose components are stage-specifically regulated at post-translational levels.

Analyzing stage-enriched, functional networks in further detail using the STRING applications of the Cytoscape software package revealed a series of potential regulatory interactions (Figure S11, Table 12): A first group of protein kinase substrates is represented by protein kinases themselves with more than 30 and 100 members of this large protein family identified in the ama and pro phosphoproteomes, respectively (Table 11, sheets D and H, see Figure 5D for phosphorylation pattern of some members of the MAP kinase pathway). An additional group of protein kinase substrates in both amastigote and promastigote phosphoproteomes include stage-specific phosphatases (18 in ama and 55 in pro) that can counteract protein kinase activities and even deactivate signaling cascades by dephosphorylation of protein kinase substrates (Figure S11, Table 11, sheets D and H). Furthermore, the stage-specific phosphoproteomes showed enrichment in functional networks linked to ribosomal biogenesis and RNA processing suggesting that protein kinases may control the establishment of stage-specific ribosomes. Finally, the identification of networks linked to proteasomal protein degradation in both ama and pro datasets not only further sustains the role of protein turnover in Leishmania stage differentiation, but uncovers a possible feedback loop between proteasomal activities controlling stage-specific degradation of protein kinases (see Figure 4), which in turn may control stage-specific ubiquitin ligases and deubiquitinating enzymes via differential phosphorylation (see Figure 5D for phosphorylation pattern of selected members of the Ubiquitin proteasomal system).

Discussion

By applying a 5-layer systems analysis to spleen-derived bona fide amastigotes and culture-derived promastigotes, we identified a series of functionally related gene clusters that are co-regulated during Leishmania donovani stage differentiation. This co-regulation occurs at multiple levels, including transcript stability, proteasomal protein turnover, and protein phosphorylation. Our study underscores the complex molecular architecture of the Leishmania developmental process and provides novel insights into possible regulatory feedback loops that may coordinate stage-specific transitions. These results raise important questions regarding the integration of these regulatory networks during differentiation, and their contributions to environmental sensing and adaptive evolution. In the following we discuss these putative networks in the context of the current literature and propose experimental strategies for their downstream, mechanistic analyses.

A first series of possible regulatory networks proposed by our results could rely on recursive or self-controlling interactions, where components of a given pathway are regulated by the pathway itself (Figure 6A). In promastigotes for example, a co-regulated, functional gene cluster showing stage-specific increase in transcript abundance was defined by the GO term ‘post-transcriptional regulation of gene expression’. In the absence of stage-specific gene dosage changes (Figure 1) and the lack of transcriptional regulation of gene expression in Leishmania (104) the stage-specific expression changes of these transcripts thus is likely regulated itself at the post-transcriptional level. Such auto-regulatory feedback is common for RNA-binding proteins (RBPs), which play a central role in controlling mRNA stability, often targeting their own transcript (105). This helps keep protein levels stable (negative feedback) or create on/off switches for developmental processes (positive feedback) (106). In T. brucei for example, the stage-specific expression of RBPs, such as ZC3H22 or RBP9, is regulated at post- transcriptional levels (107), while other RBPs, such as TbZFP3 or RBP10, can control developmental transitions by stabilizing cohorts of mRNAs, thus defining developmental regulons (108, 109). In Leishmania, RBPs are likely to participate in comparable regulatory networks. These can be investigated through methods such as crosslinking and immunoprecipitation (CLIP) to map RBP binding to 3’ untranslated regions (3’UTRs), computational approaches to identify conserved regulatory sequence motifs, and loss-of-function analyses targeting either trans-acting factors (RBPs and their interaction partners) or cis-acting elements (e.g., RBP binding sites within 3’UTRs).

Figure 6:

Our data propose a similar recursive feedback loop that seems to act on the level of protein translation, which was revealed by a paradox: while amastigotes showed increased abundance for 38 transcripts annotated for the GO term ‘ribosome biogenesis’ (Table 4, sheet E), the corresponding proteins showed reduced abundance in amastigotes compared with promastigotes (Figure 3C). This dissociation between RNA and protein abundances indicates that ribosomal components themselves may be regulated at translational levels. Increased ribosomal protein abundance in promastigotes could satisfy the requirement of increased translation capacity at this fast-growing stage (31, 110). In contrast, increased mRNA abundance of ribosomal components in amastigotes may indicate stage-specific mechanisms of mRNA storage (111) that may allow to jump start promastigote differentiation and accelerated growth following exposure to low temperature and neutral pH. This possibility is supported by our proteomics investigations of the early differentiation process, revealing ribosomal biogenesis as one of the first pathways induced over the initial 18h during the amastigote-to-promastigote developmental transition (Figure S6).

Differentiation to promastigotes was associated not only with quantitative changes of ribosomal components as indicated by co-regulated expression of ribosomal proteins and translation factors (see Table 8), but also seems to affect ribosome structure itself as indicated by stage-specific expression changes of gene sets enriched in the GO terms ‘RNA methylation’ (8 methyltransferases, including the rRNA methyltransferases Ld1S_350672300 and Ld1S_350694000) and ‘pseudouridine synthesis’ (6 pseudouridine synthases and the H/ACA snoRNP Nop10 (Ld1S_360722900) (Table 8, sheets D to H). The differential expression of snoRNAs (see Figure 2C and Table 5) and the developmentally regulated changes in rRNA modification that can affect ribosome activity (34) (Figure 2D) indeed support the existence of such stage-specific ribosomes. These results are in line with our recent reports providing indication for specialized ribosomes in Leishmania adaptation, in which we correlated changes in rRNA pseudouridylation to L. donovani fitness gain in culture (31, 66). Cryo-EM analysis of purified ribosomes from amastigotes and derived promastigotes, combined with targeted deletion of individual snoRNA genes, will be critical to uncover stage-specific structural and functional adaptations of such stage-specific ribosomes.

Finally, a third possible recursive feedback loop seems to act at the post-translational level involving members of various protein kinase (PK) families, which not only showed stage-specific expression, but also represented a main class of phosphorylation substrates themselves. This interaction may reflect stage-specific signaling cascades, where downstream kinases are activated through phosphorylation by upstream kinases, such as exemplified by the MAP kinase pathway (103, 112). Indeed, 11 members of the MAP kinase family showed stage-specific changes in their phosphorylation pattern (Figure 5D), 9 of which were only detected in the promastigote phospho-proteome and were previously associated with this stage, being implicated in flagellar biogenesis (MPK9 and MPK3) or metacyclogenesis (MPK4) (98, 102, 113). In contrast, the two MAP kinases MPK2 and MPK10 - previously linked to amastigote differentiation and virulence (99, 102) - showed an amastigote-specific increase in protein abundance and relative phosphorylation change for defined phosphorylation sites. In L. mexicana, such a recursive PK feedback loop has been identified between the interacting protein kinases LmxMKK and LmxMPK3, which cross-phosphorylate each other to regulate flagellar length (113). To further dissociate the complex regulatory relationships between parasite PKs, protein-protein interaction maps should be established (e.g. co-IP, proximity labelling), combined with in vitro kinase assays to reveal direct kinase-substrate relationships.

Most PKs have pleiotropic effects across biological systems, due to their ability to phosphorylate a wide range of functionally unrelated substrates (114). Our analysis of relative phosphorylation changes normalized to protein abundance confirms such pleotropic function in Leishmania and suggests parasite PKs at the center of a second form of possible feedback loop defined by reciprocal interactions, where components of two or more regulatory pathways may cross-control each other (Figure 6A). For example, as judged by changes in phosphorylation abundance, it seems that stage-specific PK activities act on a series of phosphatases, which in return may modify kinase activities by dephosphorylation (Table 11, sheets D and G). Such reciprocal regulation between both enzyme families is well-established across many eukaryotic systems, for example between mammalian ERK1/2 and the MAPK phosphatase DUSP6/MKP-3 (115), or yeast CDK1 and the phosphatase Cdc25, controlling mitotic entry (116). In addition, 47 differentially phosphorylated proteins were associated with ‘Translation’, some of which have been linked in Leishmania to translational control, cell cycle regulation or infectivity (29, 117–121). This opens the possibility of yet another reciprocal feedback loop, where phosphorylation of ribosomal components could fine tune translational control, which in turn may affect the expression of protein kinases themselves (122–130).

Finally, our results suggest the possibility of a third reciprocal feedback loop between protein kinases and proteins of the proteasomal system that may govern Leishmania stage transitions. Previous genetic analyses demonstrated an essential role for various components of this system in promastigote-to-amastigote transition and intracellular survival, including L. mexicana ubiquitin conjugating (E2) enzymes (UBC1/CDC34, UBC2 and UEV1), the E3 ubiquitin ligase HECT2 or various deubiquitinases (e.g. DUBs 4, 7, 13) (91, 92). By applying the irreversible proteasome inhibitor lactacystin, we extended this role to the reverse process - amastigote-to-promastigote differentiation - demonstrating that proteasomal activities are critical for both directions of the developmental cycle (Figure 4). The constitutive expression of many proteasomal components across both amastigote and promastigote stages supports a role of phosphorylation in regulating stage-specific proteasomal activities, considering that phosphorylation of DUBs has been linked to changes in their stability, localization, specificity or catabolic activity (131, 132). We indeed identified 11 DUBs and 10 ubiquitin ligases/transferases showing stage-specific phosphorylation at specific residues (Table 11, sheets D and I). This correlated with the stage-specific proteasomal degradation of 14 PKs (Figure 4F, Table 9, sheet G), including two previously implicated in Leishmania in vivo fitness (Repressor of differentiation kinase 1, Ld1S_320540600; differentiation inhibitory kinase, Ld1S_280397500) (102). These findings support the existence of a possible kinase/proteasome reciprocal feedback loop in early parasite development. The experimental validation of such regulatory interactions will require in-depth biochemical, molecular, and systems biology approaches involving global kinase-substrate mapping (114), in vitro reconstitution assays using recombinant proteins to demonstrate reciprocal kinase and proteolytic activities, and genetic perturbation test to reveal the biological impact of the regulatory interaction.

The complexity of regulatory interactions between protein kinases and the proteasomal pathway is further increased by observations showing that substrate phosphorylation can mediate subsequent ubiquitination and degradation, while ubiquitination of MAP kinase family members can regulate kinase activity and localization rather than promote degradation (139–141). Future experiments simultaneously analyzing changes in both phosphoproteome and ubiquitome during Leishmania stage-differentiation and making precision mutants that lack the ubiquitination or phosphorylation sites may help elucidate this complex regulatory interaction and identify individual proteins and their PTMs that drive parasite development

In conclusion, our data represent an important novel resource to experimentally assess the regulatory landscape of Leishmania stage development. We propose genetic feedback control as a central mechanism in parasite differentiation, with stage-specific interactions between genes and their products likely adapting the parasite phenotype to the environmental changes encountered inside its vertebrate and invertebrate hosts. In contrast to the major phenotypic shift observed between mammalian-stage amastigotes and insect-stage promastigotes, each developmental stage is stably maintained within its respective host even in the wake of environmental fluctuations. Such phenotypic robustness is known to rely on highly inter-connected, stage-specific regulatory networks (142), which in the case of Leishmania could implicate recursive and reciprocal genetic interactions between protein translation, protein degradation and protein phosphorylation (Figure 6B and C). Our results define Leishmania differentiation as an excellent model system to study how alternations between phenotypic plasticity (i.e. developmental transitions) and robustness (stable maintenance of stages) shape microbial fitness, and provides a powerful experimental framework for future studies on the emergence of stable, adaptive biological traits through feedback regulation.

Sup figure legends

Figure S1:

Figure S2:

Figure S3:

Figure S4:

Figure S5:

Figure S6:

Figure S7:

Figure S8:

Figure S9:

Figure S10:

Figure S11:

Data availability

The mass spectrometry proteomics data were deposited to the ProteomeXchange Consortium via the PRIDE partner repository with the dataset identifiers PXD035697 (phosphoproteome data) and PXD035698 (proteome upon lactacystin treatment). The DNAseq data have been deposited in NCBI's Gene Expression Omnibus and are accessible in the BioProject PRJNA1231373 (https://www.ncbi.nlm.nih.gov/bioproject/PRJNA1231373). The RNAseq data have been deposited in ArrayExpress (EMBL-EBI) (69) under the accession number E-MTAB-16528. The mass spectrometry metabolomic data are available on Figshare at https://doi.org/10.6084/m9.figshare.31027849. The scripts for GIP and giptools are available at https://github.com/susefranssen/Global_genome_diversity_Ldonovani_complex.

Acknowledgements

We thank the CEA-CNRGH for its contribution to the sequencing costs and all the CEA-CNRGH staff who performed sample preparation and sequencing for their excellent technical assistance. We acknowledge the help of the HPC Core Facility of the Institut Pasteur for this work. We thank Georges Haustant and Caroline Proux from the Biomics Platform, C2RT, Institut Pasteur, Paris, France for the library preparation and the sequencing, Hugo Varet and Rachel Legendre from Institut Pasteur, Université Paris Cité, Biostatistics and Bioinformatics Hub, Paris, France for their contribution to the RNA seq data analysis. We thank Karim Sébastien from the animal facility of the Institut Pasteur for his commitment in the care of our animals.

This work was supported by the Agence Nationale pour la Recherche Labex ‘Integrative Biology of Emerging Infectious Diseases’ contract ANR-10-LABX-62-IBEID and Labex ‘French Alliance for Parasitology and Health Care’ contract ANR-11-LABX-0024 (GFS, PP), the France Génomique National infrastructure, funded as part of the « Investissements d’Avenir » program managed by the Agence Nationale pour la Recherche contract ANR-10-INBS-09 (GH, HV, RL, CP), and the ERC SYNERGY project DecoLeishRN, Grant agreement ID: 101071613.

Additional files

Supplementary data.

Table 1.

Table 2.

Table 3.

Table 4.

Table 5.

Table 6.

Table 7.

Table 8.

Table 8b.

Table 9.

Table 10.

Table 11.

Table 12.

Table 13.

Additional information

Funding

Agence Nationale de la Recherche (ANR) (ANR-10-LABX-62-IBEID)

• Pascale Pescher

Agence Nationale de la Recherche (ANR) (ANR-11-LABX-0024)

• Pascale Pescher

Agence Nationale de la Recherche (ANR) (ANR-10-INBS-09)

• Thomas Cokelaer

EC | European Research Council (ERC) (10.3030/101071613)

• Pascale Pescher