domain: info sdv

text:
Laura Piel

K Shanmugha Rajan

Giovanni Bussotti

Hugo Varet

Rachel Legendre

Caroline Proux

Thibaut Douché

Quentin Giai-Gianetto

Thibault Chaze

Barbora Vojtkova

Nadav Gordon-Bar

Tirza Doniger

Smadar Cohen-Chalamish

Praveenkumar Rengaraj

Céline Besse

Anne Boland

Jovana Sadlova

Jean-François Deleuze

Mariette Matondo

Ron Unger

Petr Volf

Shulamit Michaeli

Pascale Pescher
email: pascale.pescher@pasteur.fr

Gerald F Späth
email: gerald.spaeth@pasteur.fr

Post-transcriptional regulation of Leishmania fitness gain

Keywords: Leishmania, experimental evolution, fitness, genomics, transcriptomics, posttranscriptional adaptation, small nucleolar RNAs

Introduction

Parasitic protozoa of the genus Leishmania are the etiologic agents of a spectrum of severe diseases known as leishmaniases that cause substantial human morbidity and are among the five most serious parasitic diseases worldwide. Today, almost 1 billion people are at risk of Leishmania infection in close to 100 endemic countries throughout tropical and subtropical regions, with over 12 million people diagnosed with the infection [START_REF] Who | Leishmaniasis in high-burden countries: an epidemiological update based on data reported in 2014[END_REF]. Leishmaniasis represents a global public health challenge: recurrent epidemics are observed in South America, the Maghreb, Middle East, East Africa and India, and Leishmania infection has been declared an emerging disease in the EU and South East Asia [START_REF] Who | Leishmaniasis in high-burden countries: an epidemiological update based on data reported in 2014[END_REF][START_REF] Lindgren | Public health. Monitoring EU emerging infectious disease risk due to climate change[END_REF].

In absence of treatment, visceral leishmaniasis (VL -also known as Kala Azar) is the most severe and fatal form of the disease, caused either by Leishmania (L.) donovani or L. infantum.

Most Leishmania species show a digenic life cycle comprising two major developmental stages that infect two distinct hosts. The motile, extracellular promastigote form of Leishmania proliferates inside the digestive tract of the sand fly insect vector, while the non-motile, intracellular amastigote form develops and proliferates inside fully acidified, macrophage phagolysosomes of mammalian hosts. Aside stage differentiation, the success of Leishmania as a pathogenic microbe relies on its capacity to adapt to a variety of environmental fluctuations encountered in their hosts via an evolutionary process. Evolutionary adaptation relies on the classical Darwinian paradigm, where spontaneous mutations and stochastic changes in gene expression generate genetically and phenotypically heterogenous populations that compete for reproductive success in a given environment, thus driving natural selection of the fittest individuals [START_REF] Koonin | Evolution of microbes and viruses: a paradigm shift in evolutionary biology?[END_REF]. While this process is well understood in viral and bacterial infections, only little information is available on evolutionary adaptation of eukaryotic pathogens, notably protozoan parasites. This is especially relevant to trypanosomatids, which -in contrast to classical eukaryotes -do not regulate expression of protein coding genes by transcriptional control. Transcription of protein coding genes in these earlybranching eukaryotes is constitutive, with genes being arranged in long, polycistronic transcription units, and mature mRNAs being generated from precursors via a highly parasite-specific, transsplicing process that is linked to polyadenylation [START_REF] Michaeli | Trans-splicing in trypanosomes: machinery and its impact on the parasite transcriptome[END_REF][START_REF] Clayton | Gene expression in Kinetoplastids[END_REF]. In the absence of classical transcriptional regulation [START_REF] Clayton | Gene expression in Kinetoplastids[END_REF], Leishmania has evolved and emphasized other forms of control, notably regulation of RNA abundance by post-transcriptional regulation and gene dosage variation [START_REF] Sterkers | The Leishmania chromosome lottery[END_REF][START_REF] Rogers | Chromosome and gene copy number variation allow major structural change between species and strains of Leishmania[END_REF][START_REF] Dumetz | Modulation of Aneuploidy in Leishmania donovani during Adaptation to Different In Vitro and In Vivo Environments and Its Impact on Gene Expression[END_REF][START_REF] Iantorno | Gene Expression in Leishmania Is Regulated Predominantly by Gene Dosage[END_REF][START_REF] Prieto Barja | Haplotype selection as an adaptive mechanism in the protozoan pathogen Leishmania donovani[END_REF].

A hallmark of Leishmania biology is the intrinsic plasticity of its genome, with frequent copy number variations (CNVs) of individual genes or chromosomes linked to drug resistance or changes in tissue tropism [START_REF] Rogers | Chromosome and gene copy number variation allow major structural change between species and strains of Leishmania[END_REF][START_REF] Brotherton | Proteomic and Genomic Analyses of Antimony Resistant Leishmania infantum Mutant[END_REF][START_REF] Downing | Whole genome sequencing of multiple Leishmania donovani clinical isolates provides insights into population structure and mechanisms of drug resistance[END_REF][START_REF] Laffitte | Plasticity of the Leishmania genome leading to gene copy number variations and drug resistance[END_REF][START_REF] Leprohon | Gene expression modulation is associated with gene amplification, supernumerary chromosomes and chromosome loss in antimony-resistant Leishmania infantum[END_REF][START_REF] Ubeda | Genome-wide stochastic adaptive DNA amplification at direct and inverted DNA repeats in the parasite Leishmania[END_REF][START_REF] Zhang | Genetic analysis of Leishmania donovani tropism using a naturally attenuated cutaneous strain[END_REF]. Combining experimental evolution and comparative genomics approaches, we recently linked both forms of genome instability with fitness gain in vitro. DNA read depth analysis of the genomes of L. donovani parasites evolving to higher fitness in culture identified amplification of a series of chromosomes as highly reproducible drivers of early culture adaptation [START_REF] Prieto Barja | Haplotype selection as an adaptive mechanism in the protozoan pathogen Leishmania donovani[END_REF]. Long-term adaptation in contrast was linked to the positive selection of gene copy number variants, which were amplified as part of functionally related, epistatic networks that allowed the emergence of phenotypes linked to ribosomal biogenesis, translation and proliferation (Bussotti et al., in preparation). Leishmania genomic adaptation thus occurs through a two-stage process reminiscent to other fast-growing eukaryotic cells (e.g. fungi and cancer cells [START_REF] Yona | Chromosomal duplication is a transient evolutionary solution to stress[END_REF][START_REF] Negrini | Genomic instability--an evolving hallmark of cancer[END_REF]), involving short-term adaptation by karyotypic changes and long-term adaptation through slower gene CNVs [START_REF] Prieto Barja | Haplotype selection as an adaptive mechanism in the protozoan pathogen Leishmania donovani[END_REF][START_REF] Bussotti | Leishmania Genome Dynamics during Environmental Adaptation Reveal Strain-Specific Differences in Gene Copy Number Variation, Karyotype Instability, and Telomeric Amplification[END_REF] ( Bussotti et al., in preparation).

Together these reports draw a complex picture of Leishmania fitness gain and raise a series of important new questions on (i) the nature of the genes that drive positive selection of the observed karyotypic changes during fitness gain in vitro, (ii) the potential fitness costs associated with karyotypic adaptation, and (iii) the mechanisms evolved by the parasite to harness genome instability for fitness gain and compensate for deleterious gene dosage effects. Here we combined experimental evolution and integrative systems approaches to address these questions and gain novel insight into regulatory mechanisms underlying Leishmania fitness gain during adaptation to culture. Our analyses reveal mechanisms at gene, transcript and protein levels that harness genome instability for fitness gain through gene dosage-dependent changes that affect post-transcriptional regulation and gene dosage-independent changes in epitranscriptomic control and ribosomal biogenesis.

Material and Methods

Animals and ethics statement

Six to eight-week-old, female mice (Mus musculus, C57BL/6JRj) and 5 female Golden Syrian hamsters (Mesocricetus auratus RjHan:AURA, weighting between 60 -70 g) were purchased from Janvier Laboratories. All animals were handled under specific, pathogen-free conditions in biohazard level 3

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. Amastigotes were recovered from infected hamster spleen and differentiated into promastigotes in M199 complete medium (M199, 10% FBS, 20 mM HEPES; 100 µM adenine, 2 mM L-glutamine, 10 µg/ml folic acid, 13.7 µM hemin, 4.2 mM NaHCO 3 , 1xRPMI1640 vitamins, 8 µM 6-biopterin, 100 units penicillin and 100 µg/ml streptomycin, pH 7.4) at 26°C. Promastigotes, derived from splenic amastigotes, were serially passaged once stationary phase was reached for less than 5 passages (EP, early passage) or 20 passages (LP, late passage) corresponding to ~ 20 and 190 generations, respectively. Strains issued from independent experimental evolution assays are identified by number (i.e. EP.1 and LP.1 are the strains resulting from experiment 1).

Parasite growth and determination of the generation time

Promastigotes in exponential growth phase were seeded at 2x10 6 (EP) or 1x10 5 (LP) parasites per ml in M199 complete medium. Parasites were counted every 24 hours during a three-day period at logarithmic growth phase, and the generation time was calculated according to the formula

𝑑𝑑𝑑𝑑𝑑𝑑𝑑𝑑𝑑𝑑𝑑𝑑𝑑𝑑𝑑𝑑 𝑡𝑡𝑑𝑑𝑡𝑡𝑡𝑡 = 𝑡𝑡×log( 2) log 𝐶𝐶(𝑡𝑡)-log 𝐶𝐶(𝑖𝑖)
. Experiments were performed in triplicates and statistical significance was assessed by t-test.

Ficoll gradient centrifugation for metacyclic promastigote enrichment

Promastigote cultures at day 3 of stationary culture phase were collected and adjusted to 3x10 8 cells/ml. Ficoll PM400 (GE Healthcare) was used to prepare a 20% stock solution in PBS and diluted for preparation of 10% and 5% Ficoll solutions. Four ml of 10% Ficoll were overlaid by 4 ml of 5% Ficoll and 4 ml of parasite suspension were layered on top of the Ficoll cushion. Tubes were centrifuged at 1,300 x g for 15 min at room temperature without brake. The metacyclic-enriched fractions were recovered at the interface between the 10% and 5% Ficoll layers. Parasites were washed with PBS and adjusted to the final concentration required for a given experiment.

DNA and RNA extraction and sequencing

The different growth kinetics between EP and LP parasites were considered when recovering parasites from logarithmic growth phase (i.e. using the same density), and stationary growth phase (i.e. using parasites kept for the same time at this phase) for DNA or RNA preparation. DNA was prepared from parasites in exponential culture phase. EP.1 and LP.1, LP.6 and LP.7 promastigotes were centrifuged at 1,600 x g for 10 min at room temperature. Approximately 1x10 8 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 the DNAs was controlled determining the DNA Integrity Number (DIN) analyzing 20 ng of DNA on a TapeStation 4200 (Agilent). One µg 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.

Total RNA was extracted from EP.1 and LP.1 cultures at logarithmic and stationary (3 days) culture phases or from EP.1 metacyclic-enriched parasites (EP.1 meta). Promastigotes were centrifuged at 3,000 x g for 10 min at room temperature 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 controlled using the Agilent Bioanalyzer. DNAse-treated RNA extracts were used for library preparation using the Truseq Stranded mRNA sample preparation kit (Illumina, San Diego, California) according to the manufacturer's instructions. An initial poly (A+) RNA isolation step (included in the Illumina protocol) was performed on total RNA to remove ribosomal RNA. Fragmentation was performed on the enriched fraction by divalent ions at high temperature. The fragmented RNA samples were randomly primed for reverse transcription followed by second-strand synthesis to create double-stranded cDNA fragments. No end repair step was necessary. An adenine was added to the 3'-end and specific Illumina adapters were ligated. Ligation products were submitted to PCR amplification. The obtained oriented libraries were controlled by Bioanalyzer DNA1000 Chips (Agilent, # 5067-1504) and quantified by spectrofluorimetry (Quant-iT™ High-Sensitivity DNA Assay Kit, #Q33120, Invitrogen). Sequencing was performed on the Illumina Hiseq2500 platform at the Biomics Center (Institut Pasteur, Paris, France) to generate single-ended, 130-bp reads bearing strand specificity.

For transcriptome-wide mapping of pseudouridine sites (Ψ-seq), total RNA was either untreated or treated with N-cyclohexyl-N-β-(4-methylmorpholinium) (CMC) in bicine buffer (0.17 M CMC in 50 mM bicine, pH 8.3, 4 mM EDTA, 7 M urea) at 37 o C for 20 min. Excess CMC was removed by ethanol precipitation. To remove CMC groups attached to G and U, the CMC-treated RNA was subjected to alkaline hydrolysis with Na 2 CO 3 (50 mM, pH 10.4) at 37 o C for 4h, as previously described [START_REF] Bakin | Mapping of the 13 pseudouridine residues in Saccharomyces cerevisiae small subunit ribosomal RNA to nucleotide resolution[END_REF][START_REF] Chikne | A pseudouridylation switch in rRNA is implicated in ribosome function during the life cycle of Trypanosoma brucei[END_REF][START_REF] Rajan | Pseudouridines on Trypanosoma brucei spliceosomal small nuclear RNAs and their implication for RNA and protein interactions[END_REF][START_REF] Liang | The spliced leader-associated RNA is a trypanosomespecific sn(o) RNA that has the potential to guide pseudouridine formation on the SL RNA[END_REF]. The reacted RNA was recovered by phenol chloroform extraction, and ethanol precipitation. An adaptor was ligated to the 3' end of the total RNA (upon fragmentation) before and after CMC treatment, and cDNA was prepared using AffinityScript reverse transcriptase (Agilent). The cDNA was then ligated to an adaptor, PCR amplified, and the samples were sequenced in an Illumina NextSeq machine in paired-end mode, 42-bp reads (20 million reads for each sample).

For the preparation of the small RNome, whole cell extracts were prepared from L. donovani EP and LP parasites (5x10 9 ), that were washed with PBS and resuspended in 20 mM Tris-HCl (pH 7.7), 25 mM KCl, and 10 mM MgCl 2, were equilibrated in a nitrogen cavitation bomb (Parr Instruments Co.) with 750 psi N2 for 1h at 4°C, and disrupted by release from the bomb. After nitrogen cavitation, the ribonucleoproteins (RNPs) were extracted with 0.4 M KCl. Ribosomes were removed by centrifugation for 3h at 35,000 rpm in a Beckman 70.1Ti rotor (150,000 x g) and the supernatant was defined as post-ribosomal supernatant (PRS). RNA was extracted after treatment with 100 µg/ml of Proteinase K, 1% SDS in the presence of 100 µg/ml DNaseI and was used for library preparation as described previously [START_REF] Chikne | A pseudouridylation switch in rRNA is implicated in ribosome function during the life cycle of Trypanosoma brucei[END_REF]. The samples were sequenced in an Illumina NextSeq machine in pairedend mode, 42-bp reads (40 million reads for each sample).

Protein extraction, digestion and LC-MS/MS acquisition

Exponentially growing EP.2, EP.3, EP.4, EP.5 and LP.2, LP.3, LP.4 and LP.5 promastigotes were centrifuged at 1600 x g for 10 min at 4°C and washed three times with cold PBS. Parasite lysates were prepared in 8 M urea, 50 mM Tris, supplemented with a protease inhibitor cocktail (cOmplete™ from Roche) and a phosphatase inhibitor cocktail (PhosStop from Roche), 1 ml of lysis buffer per 1.5x10 9 promastigotes. After 10 min incubation at 4°C followed by sonication for 5 min (sequence of 10s pulse and 20s pause) the lysates were centrifuged 15 min at 14,000 x g, 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 a control of the protein pattern of all the extracts was performed by SDS-PAGE and silver staining. All the biological samples were further processed for MS-based proteomics approach, data acquisition, and statistical analyses.

Biological samples were adjusted to 1.3 µg.µl -1 in lysis buffer. Disulfide bridges were reduced in 5 mM DTT (Sigma -43815) for 30 min and alkylated in 20 mM iodoacetamide (Sigma -I1149) for 30 min at room temperature in the dark. Protein samples were diluted 10 fold in 50 mM Tris-HCl and digested with Sequencing Grade Modified Trypsin (Promega -V5111) at a Protein:Trypsin ratio 50:1 overnight. Then a second digestion was performed to complete this step. Proteolysis was stopped by adding formic acid (FA, Fluka -94318) at a 1% final concentration. Resulting peptides were desalted using Sep-Pak SPE cartridge (Waters) according to the manufacturer's instructions. Peptides were concentrated to almost dryness and were resuspended in 2% Acetonitrile (ACN) / 0.1% FA just before LC-MS/MS injection.

All analyses were performed on a Q Exactive TM Plus Mass Spectrometer (Thermo Fisher Scientific) coupled with a Proxeon EASY-nLC 1000 (Thermo Fisher Scientific). One µg of peptides was injected into a home-made 50 cm C18 column (1.9 μm particles, 100 Å pore size, ReproSil-Pur Basic C18, Dr.

Maisch GmbH, Ammerbuch-Entringen, Germany). Column equilibration and peptide loading were performed at 900 bars in buffer A (0.1% FA). Peptides were separated with a multi-step gradient of 2 to 5% buffer B (80% ACN, 0.1% FA) for 5 min, 5 to 22% buffer B for 150 min, 22 to 45% buffer B for 60 min, 45 to 80% buffer B for 10 min at a flow rate of 250 nL/min over 240 min. Column temperature was set to 60°C. MS data were acquired using Xcalibur software using a data-dependent method. 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 3E 6 , 20 ms and 1E 6 , 60ms respectively. An automatic selection of the 10 most intense precursor ions was activated (Top 10) with a 45s dynamic exclusion. The isolation window was set to 1.6 m/z and normalized collision energy fixed to 28 for HCD fragmentation. We used an underfill ratio of 1.0% corresponding to an intensity threshold of 1.7E 5 . Unassigned precursor ion charge states as well as 1, 7, 8 and >8 charged states were rejected and peptide match was disabled.

Data analyses

WGS analysis: Genomic DNA reads were aligned to the L. donovani Ld1S reference genome (Bussotti & Myler unpublished, https://www.ncbi.nlm.nih.gov/bioproject/PRJNA396645) 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 [START_REF] Li | The Sequence Alignment/Map format and SAMtools[END_REF].

RealignerTargetCreator and IndelRealigner from the GATK suite were run to homogenize indels [START_REF] Depristo | A framework for variation discovery and genotyping using nextgeneration DNA sequencing data[END_REF].

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 [START_REF] Quinlan | BEDTools: a flexible suite of utilities for comparing genomic features[END_REF]. 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." The chromosome sequencing coverage was used to evaluate aneuploidy between EP.1 and LP.1 samples. For this purpose the nucleotide coverage was normalized by the median genomic coverage. Then for each sample and for each chromosome, the median sequencing coverage was computed for contiguous windows of 2,500 bases. Eventually the window coverage was further normalized by chromosome 36 median coverage and multiplied by two. The gene counts were produced using featureCounts (version 1.4.6-p3 [START_REF] Liao | featureCounts: an efficient general purpose program for assigning sequence reads to genomic features[END_REF]) with these parameters: -s 0 -t gene -g gene_id and were normalized according to the median-ratio method.

RNAseq analysis: For total RNAseq data, the bioinformatics analysis was performed using the RNAseq pipeline from Sequana [START_REF] Cokelaer | Sequana': a Set of Snakemake NGS pipelines[END_REF]. Reads were cleaned of adapter sequences and low-quality sequences using cutadapt version 1.11 [START_REF] Martin | Cutadapt removes adapter sequences from high-throughput sequencing reads[END_REF]. Only sequences of at least 25 nucleotides in length were considered for further analysis. STAR version 2.5.0a [START_REF] Dobin | STAR: ultrafast universal RNA-seq aligner[END_REF], with default parameters, was used for alignment on the reference genome (GCA_002243465.1). Genes were counted using featureCounts version 1.4.6-p3 [START_REF] Liao | featureCounts: an efficient general purpose program for assigning sequence reads to genomic features[END_REF] from Subreads package (parameters: -t gene -g gene_id -s 1). Count data were analyzed using R version 3.6.1 [START_REF] Team | R: A Language and Environment for Statistical Computing[END_REF] and the Bioconductor package DESeq2 version 1.24.0 [START_REF] Love | Moderated estimation of fold change and dispersion for RNA-seq data with DESeq2[END_REF]. The normalization and dispersion estimation were performed with DESeq2 using the default parameters and statistical tests for differential expression were performed applying the independent filtering algorithm. A generalized linear model including the replicate effect as blocking factor was set in order to test for the differential expression between EP.1 log, EP.1 stat, EP.1 meta, LP.1 log and LP.1 stat.

For each pairwise comparison, raw p-values were adjusted for multiple testing according to the Benjamini and Hochberg (BH) procedure [START_REF] Benjamini | Controlling the False Discovery Rate: A Practical and Powerful Approach to Multiple Testing[END_REF] and genes with an adjusted p-value lower than 0.05 were considered differentially expressed. The RNAseq data have been deposited in NCBI's Gene Expression Omnibus (Edgar et al., 2002) and are accessible through GEO Series accession number GSE165615 (https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE165615).

For small RNome analysis, Ψ-seq and detection of pseudouridylated sites, the 42 bp sequence reads obtained from the Illumina Genome Analyzer were first trimmed of Illumina adapters using the FASTX toolkit (http://hannonlab.cshl.edu/fastx_toolkit), and reads of 15 nucleotides or less were discarded from subsequent analysis. The remaining reads were mapped to the reference genome (GCA_002243465.1) using SMALT v0.7.5 (https://www.sanger.ac.uk/tool/smalt-0/) with the default parameters. Only properly paired partners were retained. Each read pair was "virtually" extended to cover the area from the beginning of the first read to the end of its partner. For each base, the number of reads initializing at that location as well as the number of reads covering the position were calculated. A combination of BEDTools v2.26.0 Suite (http://bedtools.readthedocs.io/en/latest/) and in-house Perl scripts was used to calculate the Ψratio and Ψ-fc (fold change), as previously described [START_REF] Chikne | A pseudouridylation switch in rRNA is implicated in ribosome function during the life cycle of Trypanosoma brucei[END_REF][START_REF] Rajan | Pseudouridines on Trypanosoma brucei spliceosomal small nuclear RNAs and their implication for RNA and protein interactions[END_REF].

Proteomics analysis: Raw data were analyzed using MaxQuant software version 1.5.3.8 [START_REF] Cox | MaxQuant enables high peptide identification rates, individualized p.p.b.-range mass accuracies and proteome-wide protein quantification[END_REF] using the Andromeda search engine [START_REF] Cox | Andromeda: a peptide search engine integrated into the MaxQuant environment[END_REF]. The MS/MS spectra were searched against the Ld1S database (Bussotti et al., in preparation). The settings for the search included (i) trypsin digestion with a maximum of two missed cleavages, (ii) variable modifications for methionine oxidation and Nterminal acetylation, and (iii) fixed modification for cysteine carbamidomethylation. The minimum peptide length was set to 7 amino acids and the false discovery rate (FDR) for peptide and protein identification was set to 0.01. The main search peptide tolerance was set to 4.5 ppm and to 20 ppm for the MS/MS match tolerance. The setting 'second peptides' was enabled to identify cofragmentation events. Quantification was performed using the XIC-based Label-free quantification (LFQ) algorithm with the Fast LFQ mode as previously described [START_REF] Cox | Accurate proteome-wide label-free quantification by delayed normalization and maximal peptide ratio extraction, termed MaxLFQ[END_REF]. Unique and razor peptides, including modified peptides, with at least two ratio counts were accepted for quantification. The mass spectrometry proteomics data were deposited to the ProteomeXchange Consortium via the PRIDE partner repository with the dataset identifier PXD020236 [START_REF] Perez-Riverol | The PRIDE database and related tools and resources in 2019: improving support for quantification data[END_REF].

For the differential analyses, proteins categorized as 'reverse', 'contaminant' and 'only identified by site' were discarded from the list of identified proteins. After log2 transformation, LFQ values were normalized by median centering within conditions (normalizeD function of the R package DAPAR [START_REF] Wieczorek | DAPAR & ProStaR: software to perform statistical analyses in quantitative discovery proteomics[END_REF]). Remaining proteins without any LFQ value in one of the conditions (either EP or LP) and at least two values in the other condition were considered as exclusively expressed proteins. Missing values across the four biological replicates were imputed using the imp.norm function of the R package norm (norm: Analysis of multivariate normal datasets with missing values. 2013 R package version 1.0-9.5). A limma t test was applied to determine proteins with a significant difference in abundance while imposing a minimal fold change of 2 between the conditions to conclude that they are differentially abundant [START_REF] Smyth | limma: Linear Models for Microarray Data[END_REF][START_REF] Ritchie | limma powers differential expression analyses for RNA-sequencing and microarray studies[END_REF]. An adaptive Benjamini-Hochberg procedure was applied on the resulting p-values using the function adjust.p of R package cp4p [START_REF] Giai Gianetto | Calibration plot for proteomics: A graphical tool to visually check the assumptions underlying FDR control in quantitative experiments[END_REF] and the robust method described in Pounds et al. [START_REF] Pounds | Robust estimation of the false discovery rate[END_REF] to estimate the proportion of true null hypotheses among the set of statistical tests. The proteins associated to an adjusted p-value inferior to a False Discovery Rate (FDR) of 0.01 have been considered as significant and differentially abundant proteins.

Gene Ontology (GO)-enrichment analyses:

The Biological Networks Gene Ontology tool (BiNGO) plugin of the Cytoscape software package (version 3.8.2) was used. A Benjamini & Hochberg false discovery rate with a significance level of 0.05 was applied. The lists of L. donovani GO terms were built in house. In order to assign the Gene Ontology Identifiers (GO IDs) we combined the GO-derived identifiers with the ones available from the corresponding orthologs in target species: LdBPK, L. infantum, L. major, L. mexicana, Typanosoma brucei brucei 927 (Tbru) and Typanosoma cruzi (Tcru).

For each target species we retrieved both the "curated" and "computed" GO IDs from TriTrypDB on the 11/09/2019. OrthoFinder with the DIAMOND search program was applied to establish orthology between the genes in Ld1S and in target species. In "one-to-many" orthology relations we concatenated all the non-redundant GO IDs from all the homologs. The GO IDs were then assigned based on the hierarchy:

LdBPK curated > LdBPK GO > L. infantum curated > L. major curated > L. mexicana curated > Tbru curated > Tcru curated > LdBPK computed > L. infantum computed > L.
major computed > L. mexicana computed > Tbru computed > Tcru computed. The GO IDs were assigned if not present in any higher rank GO ID data set. The GO IDs of snoRNAs, UsnRNA, SLRNA and 7SL classes defined by homology with L. major Friedlin genes were manually attributed. Overall, we assigned biological process (BP), molecular function (MF) and cellular component (CC) GO IDs to 5,246, 4,521 and 7,236 Ld1S genes.

Northern blot analyses

Total RNA extracted from EP and LP cells (10 µg) were separated on 10% acrylamide denaturing gels, transferred to nitrocellulose membranes and analyzed by autoradiography. RNA probes were prepared by in vitro transcription using α-32 P-UTP [START_REF] Chikne | A pseudouridylation switch in rRNA is implicated in ribosome function during the life cycle of Trypanosoma brucei[END_REF]. Three independent northern blots were performed.

Bone Marrow Derived Macrophages and infection

Bone marrow exudate cells were recovered from tibias and femurs of C57BL/6JRj female mice (Janvier Labs) and macrophages differentiated in DMEM complete medium (DMEM, 15% FBS, 10 mM HEPES, 50 µM 2-mercaptoethanol, 50 units of penicillin and 50 µg/ml of streptomycin) supplemented with 75 ng/ml of recombinant mouse colony stimulating factor-1 (rmCSF-1, ImmunoTools) [START_REF] Aulner | High content analysis of primary macrophages hosting proliferating Leishmania amastigotes: application to anti-leishmanial drug discovery[END_REF]. A total of 1.5x10 5 Bone Marrow Derived Macrophages (BMDMs) were plated on glass coverslips in 24wells plates and incubated overnight at 37°C, 5% CO 2 prior to Leishmania infection.

Promastigotes from stationary culture phase or metacyclic-enriched parasite fractions were pelleted by centrifugation at 3,000 x g for 10 min at room temperature and re-suspended in PBS. The concentration was adjusted to 6x10 7 parasites per ml and 50 µl were added to the BMDM cultures at a multiplicity of infection (MOI) of 20 parasites per 1 macrophage. Plates were centrifuged at 300 x g for 5 min at room temperature to allow for a faster sedimentation of the parasites onto the macrophage monolayer. After 2h of contact, coverslips were washed by successive baths in prewarmed PBS to remove extracellular parasites and transferred into new 24-wells plates containing fresh pre-warmed DMEM culture medium supplemented with 30 ng/µl of rmCSF-1. At 4, 24, 48 and 168h post-infection, cells were fixed in 4% paraformaldehyde (Electron Microscopy Science) and macrophage and parasite nuclei were stained with Hoechst 33342. Images were acquired using a

Zeiss Apotome microscope at 40x magnification connected to an Axiocam camera. All the infections were performed in triplicates and at least 100 macrophages were counted per coverslip. The total numbers of infected and non-infected macrophages were recorded and the percentage of infection, the number of parasites per 100 cells and the number of parasites per infected macrophages was calculated and normalized to the values obtained at the initial 4-hour time point. The experiment was performed three times.

Morphological analyses

Parasites were seeded on poly L-lysine treated coverslips and fixed in 2.5% glutaraldehyde. Coverslips were mounted on glass slides using Mowiol® 4-88 (Sigma-Aldrich). Images were acquired using an Axiophot microscope at 63x magnification and an Andor camera. Length and width of the parasite cell body, and flagellum length were measured for at least 200 promastigotes using the Image J Fidji software package (https://imagej.net/). The ratios flagellum over body length and body length over body width were determined for the 200 parasites and the Kruskal-Wallis test was used for statistical analysis. The experiment was performed in duplicate.

Sand fly infection

The colony of Phlebotomus orientalis (originating from Ethiopia), the natural vector of L. donovani, was maintained in the insectary of the Department of Parasitology, Charles University in Prague, under standard conditions (26°C on 50% sucrose and 14h light/10h dark photoperiod) as described previously [START_REF] Volf | Establishment and maintenance of sand fly colonies[END_REF][START_REF] Seblova | Phlebotomus orientalis sand flies from two geographically distant Ethiopian localities: biology, genetic analyses and susceptibility to Leishmania donovani[END_REF].

Promastigotes from logarithmic-phase cultures (day 3-4 in culture) were washed twice in saline solution and resuspended in heat-inactivated rabbit blood at a concentration of 1x10 6 promastigotes/ml. Sand fly females (5-9 days old) were infected by feeding through a chick-skin membrane (BIOPHARM) on a promastigote-containing suspension. Engorged sand flies were separated and maintained under the same conditions as the colony. On day 8 post blood meal (PBM), 150 sand fly females were dissected. The thoracic parts and abdominal parts of infected guts were collected separately and pooled together into two samples: thoracic parts of gut (TP) and abdominal parts of gut (AP). The exact numbers of all parasite stages were calculated using a Burker apparatus and the proportion of metacyclic forms was identified on a Giemsa-stained smears separately for TP and AP. Leishmania with flagellar length < 2 times body length were scored as procyclic forms and those with flagellar length ≥2 times body length as metacyclic forms [START_REF] Sadlova | The stage-regulated HASPB and SHERP proteins are essential for differentiation of the protozoan parasite Leishmania major in its sand fly vector, Phlebotomus papatasi[END_REF].

Results

L. donovani culture adaptation causes a fitness tradeoff between proliferation and infectivity.

In microbial culture, fitness gain (defined as reproductive capacity) largely equals the level of cell proliferation. Adaptation to in vitro growth thus represents a simple experimental system to assess mechanisms underlying fitness gain. Here we applied such an experimental evolution approach on L.

donovani amastigotes isolated from infected hamster spleen. Derived promastigotes at early-passage (EP.1) and late-passage (LP.1) were monitored for growth and infectivity with the aim to assess regulatory mechanisms underlying fitness gain and fitness cost observed during culture adaptation.

Analyzing cell growth during promastigote culture adaptation revealed robust fitness gain as judged by the reduction in generation time from 13.76 +/-1.18 hours for EP.1 to 9.76 +/-0.93 hours for LP.1 promastigotes (Figure 1A). We next evaluated fitness of these parasites in intracellular macrophage infection, where reproductive success depends on parasite resistance to host cell cytolytic activities, amastigote differentiation and proliferation. BMDMs were incubated with EP.1 and LP.1 promastigotes from day-3 stationary culture (referred to as EP.1 stat and LP.1 stat) and intracellular growth was monitored microscopically for 7 days as previously described [START_REF] Pescher | Quantitative proteome profiling informs on phenotypic traits that adapt Leishmania donovani for axenic and intracellular proliferation[END_REF]. Even though the number of EP.1 stat and LP.1 stat intracellular parasites decreased over the first 24h post-infection, only EP.1 parasites recovered and established persistent infection, while the number of LP.1 parasites steadily declined during the subsequent 6 days (Figure 1B and Figure S1A). Together these results document a fitness tradeoff occurring during L. donovani culture adaptation, with accelerated in vitro growth causing virulence attenuation.

We then tested if the fitness cost of LP.1 stat in infectivity was due to a differentiation defect of infectious metacyclic promastigotes. Ficoll gradient centrifugation [START_REF] Spath | A lipophosphoglycan-independent method for isolation of infective Leishmania metacyclic promastigotes by density gradient centrifugation[END_REF] indeed revealed a 5.5-fold reduction in the number of metacyclic parasites from 3.8% in EP.1 stat to 0.69% in LP.1 stat cultures (Figure 1C), which were in addition compromised to establish macrophage infection (Figure 1D and S1B). These results document that the fitness cost in LP.1 meta not only affects the quantity but also the quality of differentiating metacyclic parasites. This was further confirmed by their atypical morphology that was significantly different to bona fide, sand fly-isolated metacyclic promastigotes (Figure S2A), corresponding to leptomonad-like forms as judged by flagellum/body-length ratio and body shape [START_REF] Rogers | The role of promastigote secretory gel in the origin and transmission of the infective stage of Leishmania mexicana by the sandfly Lutzomyia longipalpis[END_REF][START_REF] Lei | Population changes in Leishmania chagasi promastigote developmental stages due to serial passage[END_REF] (Figure 1E, Figure S2A-C). Surprisingly, this defect in metacyclogenesis was compensated in cultures that were passaged using enriched EP.1 meta rather than EP.1 stat (Figure 1F), which likely re-enter cell cycle faster thus outcompeting metacyclic-competent parasites.

Transcriptome profiling informs on mechanisms underlying fitness tradeoff. We performed RNAseq analyses using poly(A+)-enriched mRNA obtained from three replicates of EP. 2C, left panel), with functional enrichment observed for the GO terms 'ribosome biogenesis', 'ribosome assembly', and 'rRNA processing' (Figure 2C, middle panel and Supplementary Table 2g). LP.1 log fitness gain thus likely reflects an increase in translation efficiency, which may allow for accelerated growth observed in these cells. Further interrogation combining GO analysis and manual inspection of gene annotation revealed increased abundance of other transcripts implicated in various regulatory processes linked to proliferation (Figure 2C, right panel), including epigenetic or epi-transcriptomic regulation (10 genes, e.g. Ld1S_110036500 encoding for a Pseudouridylate synthase 10, Ld1S_260334600 encoding for a RNA pseudouridylate synthase and Ld1S_330597500 encoding for a Histone methyltransferase DOT1) and cell cycle control (12 genes, e.g. Ld1S_050817000 encoding for CYC2-like cyclin, or Ld1S_330603400 encoding for the cell division control protein CDC45) (see Supplementary Table 2e for more examples).

In contrast, no GO enrichment was observed for the 437 transcripts showing reduced abundance in LP.1 log. Manual inspection of gene annotations identified various pathways implicated in metabolism and energy production (e.g. genes encoding for respiratory chain proteins, amino acid and sugar metabolism, fatty acid biosynthesis), signaling (numerous kinases and phosphatases) and flagellar biogenesis (including four genes encoding for paraflagellar rod components) (Figure S3C).

These pathways suggest a potential retooling of the LP.1 log energy metabolism in the nutrition-rich culture environment, and selection against motility, which may liberate the energy required for faster growth. Surprisingly, one of the most significant decreases in transcript abundance in these cells was observed for a gene encoding for a 5S ribosomal protein, along four other genes encoding for ribosomal components (Supplementary Table 2e), even though other ribosomal components were upregulated in LP.1 log. This result provided a first indication that LP.1 log fitness gain not only depends on the quantity, but likely also the quality or type of ribosomes, which may be adapted either for higher translation efficiency or translational control of the fitness-adapted expression profile.

Finally, we assessed changes in transcript abundance observed at stationary growth phase in LP.1 compared to EP.1 promastigotes to gain further insight into mechanisms of fitness loss (i.e. attenuated infectivity). We identified 662 transcripts with increased abundance in LP.1 stat (Figure 2D, left panel). Enrichment was observed for the GO terms 'ribosomal large subunit assembly', 'rRNA processing' and 'regulation of gene expression' (Figure 2D, middle panel). In contrast, LP.1 stat promastigotes showed reduced abundance for 710 transcripts, including several transcripts linked to the GO terms 'histone phosphorylation', 'centriole replication', 'telomere maintenance', 'reciprocal meiotic recombination' and 'response to DNA damage stimulus' (Figure 2D, right panel). Likewise, decreased abundance was observed for transcripts associated with the GO terms 'evasion or tolerance of immune response of other organism involved in symbiotic interaction', 'cellular response to osmotic stress' and 'arginine biosynthetic process' (Figure 2D, right panel). The depletion of LP.1 stat parasites in these functions could be associated with attenuated infectivity observed in LP.1 meta (see Figure 1D), which may be further aggravated by reduced expression in these parasites of 31 amastin surface glycoproteins and amastin-like proteins as well as 32 other genes previously associated with parasites infectivity, including GP63 (Supplementary Table 2f) [START_REF] Mcgwire | Genetic rescue of surface metalloproteinase (gp63)-deficiency in Leishmania amazonensis variants increases their infection of macrophages at the early phase[END_REF][START_REF] Spath | Lipophosphoglycan is a virulence factor distinct from related glycoconjugates in the protozoan parasite Leishmania major[END_REF].

In conclusion, our data link increased fitness in in vitro growth of LP.1 log to a gain-offunction phenotype associated with proliferation, ribosomal biogenesis, and translation. Conversely, the reduced fitness in infectivity of LP.1 stat was associated with a loss-of-function phenotype linked to reduced expression of virulence genes.

Post-transcriptional adaptation during promastigote fitness gain. The observed changes in transcript abundance during in vitro fitness gain may be caused by increased gene dosage due to chromosomal amplification [START_REF] Dumetz | Modulation of Aneuploidy in Leishmania donovani during Adaptation to Different In Vitro and In Vivo Environments and Its Impact on Gene Expression[END_REF][START_REF] Prieto Barja | Haplotype selection as an adaptive mechanism in the protozoan pathogen Leishmania donovani[END_REF]. Indeed, comparative genomic analysis of EP.1 and LP.1 parasites revealed aneuploidy for 9 chromosomes during culture adaptation, including trisomies for chromosomes (chr) 5, 23, 26, and 33 (Figure S4A and Supplementary Table 3a). Plotting normalized genomic versus transcriptomic read depth ratios for EP.1 and LP.1 log and stat parasites correlated 75% of the up-regulated genes in LP.1 log and 42% in the LP.1 stat promastigotes with amplified chromosomes (Figure 3A), affecting various biological processes associated with the LP.1 fitness tradeoff (Figure S5). Nevertheless, interrogating more specifically the read-depth ratios for trisomic chr 5 and 26 uncovered surprisingly high, gene dosage-independent fluctuations of RNA abundance in EP.1 and LP.1 promastigotes (Figure 3B and Supplementary Table 4d-e). While a significant fraction of transcripts on the trisomic chromosomes showed the expected 1.5-fold increase in abundance, numerous transcripts either exceeded this increase or on the contrary were expressed at lower-thanexpected abundance. Such fluctuations were also observed for the LP.1/EP.1 ratios of disomic chromosomes (see chr 36, Figure 3B as an example and Supplementary Table 4f).

We Post-transcriptional increase in mRNA abundance was observed in LP.1 log parasites for the biological processes 'rRNA processing', 'ribosome biogenesis', 'translational initiation', and 'nuclear transport' (Figure 3D and Supplementary Table 4g). In contrast, manual inspection revealed posttranscriptional decrease in abundance of mRNAs involved in flagellar biogenesis or EP.1-specific, ribosomal components (Figure 3E and Supplementary Table 4b). Significantly, reduction of both DNA and RNA read depth was observed for two NIMA-related protein kinases on chr 36 that we previously associated with in vitro fitness gain (Bussotti et al., in preparation) (see Figure 3E and Supplementary Table 4b), firmly linking the depletion in RNA and DNA abundance and eventually the deletion of these genes to accelerated growth.

In conclusion, the global analysis of the EP.1 and LP.1 transcriptomes uncovers posttranscriptional regulation as an important novel processes that may affect Leishmania fitness gain, which can likely buffer against deleterious effects of genome instability and adapt mRNA abundance in a gene dosage-independent manner to a given environment.

The fitness-adapted proteome is highly robust and enriched for GO terms associated with ribosomal biogenesis and post-transcriptional regulation. We applied a label-free, quantitative proteomics approach to assess how genomic and post-transcriptional adaptation during in vitro fitness gain impact protein abundance. Analyzing four independent, biological replicates of EP and LP strains (termed EP.2-5 and LP.2-5, Figure S6) identified a total of 6,050 proteins considering all samples, including 59 and 110 proteins that were exclusively detected in LP and EP parasites, respectively (Figure 4A and Supplementary Table 5). Considering all proteins that showed a statistically different abundance (Figure 4B), the majority of differentially expressed genes were shared in all four independent LP strains (566 of 788 total, 71%). These data reveal a surprising convergence of the fitness-adapted proteomes despite possible karyotypic variations between strains (Figure S4B), and inform on common pathways that are under convergent selection in LP strains during in vitro fitness gain. Just as observed on RNA levels, flagellar biogenesis is clearly under negative selection during culture adaptation, with reduced protein abundance observed for 38

proteins linked to motility encoded on 20 chromosomes (Figure 4C and Supplementary Table 5g).

Another key process associated with adaptation was translation: 30 proteins encoded on 11 different chromosomes were under positive selection in LP strains (e.g. various ribosomal proteins of the 39S, 40S, 60S, L22e, and S25 families, or the ribosomal assembly protein RRB1), while only two RNA binding proteins encoded on two chromosomes were under negative selection in the same parasites (Figure 4C and Supplementary Table 5c).

We next assessed the level of correlation between protein abundance, gene dosage variation and transcript abundance to gain further insight into regulatory mechanisms underlying Leishmania fitness gain. Even though the proteomics data set was obtained with four independent biological replicates (EP.2-5 and LP.2-5), the highly reproducible nature of the chr 5 and chr 26 trisomies observed in all our previous experimental evolution experiments (Figure S4B) provided a useful benchmark to assess correlations between the different data sets for at least these chromosomes.

Based on our systems comparison, three different regulatory clusters can be distinguished for chr 5 and 26: Cluster 1 includes proteins whose change in abundance correlates to gene dosage and RNA abundance (Figure 4D left and right panels, upper right and lower left quadrants), including three DNA J proteins, a HSP70 like protein and BiP, suggesting that increased stress resistance may be the driving force for the selection of these aneuploidies (Supplementary Table 6h). Cluster 2 includes proteins whose levels only correlate with mRNA abundance but not gene dosage, suggesting posttranscriptional regulation (Figure 4D right panel, see also Figure 3A and Supplementary Table 6i).

Finally, cluster 3 represents proteins whose levels do not correlate with mRNA abundance (Figure 4D left panel, lower right and upper left quadrants), which are either regulated at translational levels or by protein turn-over (Supplementary Table 6j).

Gene ontology analysis of all proteins that fall into regulatory cluster 1 (as defined by the upper right and lower left quadrants of figure 4D) revealed a strong enrichment for the term 'posttranscriptional regulation of gene expression' (Figure 4E and Supplementary Table 6b and6c). This enrichment is driven by the coordinated increase in expression of various proteins with known functions in RNA turnover (e.g. pumilio-domain protein encoded by Ld1S_330590200) and a series of proteins previously linked to post-transcriptional regulation in T. brucei such as EIF4G1 or PRP38 pre-mRNA processing factor [START_REF] Erben | A genome-wide tethering screen reveals novel potential posttranscriptional regulators in Trypanosoma brucei[END_REF][START_REF] Lueong | Gene expression regulatory networks in Trypanosoma brucei: insights into the role of the mRNA-binding proteome[END_REF][START_REF] De Pablos | The mRNA-bound Proteome of Leishmania mexicana: Novel Genetic Insight into an Ancient Parasite[END_REF]. In addition, the enrichment for the GO term 'ciliary cell motility' is driven by the under-representation of this process in the LP proteome (Figure S6G and Supplementary Table 6b and6c).

In conclusion, the Leishmania proteome undergoes reproducible, gene dosage-dependent and -independent changes during fitness gain in vitro. The robustness of proteomic adaptation indicates the presence of regulatory mechanisms that compensates for the genetic and transcriptomic variability between independent LP strains. At least under our experimental conditions, gene dosage-dependent changes modulate post-transcriptional regulation, which results in stabilization of various transcripts implicated in rRNA processing and ribosomal biogenesis. This indicates the generation of fitness-adapted ribosomes, which could generate the observed proteomic robustness. The role of ncRNAs in ribosomal biogenesis [START_REF] Rajan | Unique Aspects of rRNA Biogenesis in Trypanosomatids[END_REF] primed us in the following to carry out a dedicated small RNome analysis in EP and LP parasites to further assess the generation of specialized ribosomes.

Mapping the Leishmania non-coding transcriptome correlates snoRNA expression and rRNA modification to Leishmania fitness gain. Non-coding (nc) RNAs such as small nuclear (sn), small nucleolar (sno), ribosomal (r) or transfer (t) RNAs play essential roles in post-transcriptional regulation and translational control [START_REF] Rajan | Pseudouridines on Trypanosoma brucei spliceosomal small nuclear RNAs and their implication for RNA and protein interactions[END_REF][START_REF] Kamina | Ribosome Assembly in Trypanosomatids: A Novel Therapeutic Target[END_REF]. While our data suggested an important role of these regulatory processes in genome-independent fitness gain (see above), they did not inform on underlying ncRNAs as our RNAseq analyses used poly (A+)-enriched mRNA. We therefore performed a dedicated analysis of the small RNome in EP.1 and LP.1 L. donovani parasites. We first annotated our PacBio LD1S reference genome for ncRNAs using bioinformatics approaches (ortholog mapping, de novo annotation) as well as unmapped RNAseq reads of post-ribosomal supernatants (PRS) that are enriched in ncRNAs (Figure S7). These efforts established a very first repertoire of ncRNAs in any Leishmania species and identified 1504 genes encoding for snoRNA organized in 42 clusters on twenty four chromosomes, 83 tRNA genes, 12 snRNA genes and 140 SL RNA genes (Supplementary Table1, Figure 5A andB). Considering the trisomic chromosomes, we found 269 snoRNA genes on chr 5, 22 on chr 23, 160 on chr 26, and 193 on chr 33. We investigated more specifically the role of snoRNAs in LP.1 fitness gain given the enrichment of the fitness-adapted transcriptome in the GO terms 'ribosomal biogenesis' and 'rRNA processing' (see Figure 2C). snoRNAs guide specific modifications of rRNA, such as methylation and pseudouridylation, which in turn change the specificity of the ribosome towards certain mRNAs and thus control translation [START_REF] Watkins | The box C/D and H/ACA snoRNPs: key players in the modification, processing and the dynamic folding of ribosomal RNA[END_REF]. We prepared whole cell lysates from both EP.1 log and LP.1 log parasites, removed the ribosomes by ultracentrifugation, and prepared libraries from the post-ribosomal supernatant (PRS). From 174 detected snoRNAs, 93 showed a more than 2-fold change in LP.1 compared to EP.1, revealing a global increase of snoRNA abundance during fitness gain (Supplementary Table 7). Increased abundance was confirmed for 7 out of the 8 snoRNAs probed by Northern blot analysis of the PRS (Figure 5C andD).

Next, we examined if the increase in sno RNA abundance affected the level of rRNA pseudouridylation (Ψ) by applying a modified RNAseq protocol (termed Ψ-seq) using total RNA from EP.1 and LP.1. We detected two hyper-modified rRNA sites in all three biological replicates at positions Ψ1265 and Ψ1282 inside the peptidyl transferase center (PTC), which correlated with the increased abundance of the corresponding snoRNAs that guide these modifications (Figure 5F, Supplementary Table 7) [START_REF] Eliaz | Genome-wide analysis of small nucleolar RNAs of Leishmania major reveals a rich repertoire of RNAs involved in modification and processing of rRNA[END_REF]. Our data thus provide a first link of snoRNA expression and rRNA modification to Leishmania fitness gain, which further supports the generation of fitness-adapted ribosomes and suggests translational control -in addition to genomic and post-transcriptional adaptation -as yet another, gene dosage-independent mechanism of Leishmania evolutionary adaptation.

Discussion

A common hallmark of all microbial pathogens is their capacity to adapt to unpredictable fluctuations in their host environments through an evolutionary process, where genetically heterogenous individuals constantly compete for survival [START_REF] Lenski | What is adaptation by natural selection? Perspectives of an experimental microbiologist[END_REF]. Here we combined experimental evolution and integrative systems approaches to uncover mechanisms of fitness gain in the human pathogen Leishmania donovani. Our study provides first evidence that these parasites combine regulatory processes at genomic, post-transcriptomic and translational levels to establish highly robust fitness phenotypes while maintaining genetic heterogeneity thereby avoiding genetic death.

In the absence of classical, promoter-driven control of gene expression, Leishmania relies on alternative mechanisms to regulate transcript and protein abundance, including regulated mRNA turn over and translational control [START_REF] Clayton | Gene expression in Kinetoplastids[END_REF][START_REF] Clayton | Post-transcriptional regulation of gene expression in trypanosomes and leishmanias[END_REF]. These parasites further use a highly unusual, genomic form of gene expression regulation, where changes in chromosome and gene copy number control transcript abundance via gene dosage [START_REF] Rogers | Chromosome and gene copy number variation allow major structural change between species and strains of Leishmania[END_REF][START_REF] Dumetz | Modulation of Aneuploidy in Leishmania donovani during Adaptation to Different In Vitro and In Vivo Environments and Its Impact on Gene Expression[END_REF][START_REF] Iantorno | Gene Expression in Leishmania Is Regulated Predominantly by Gene Dosage[END_REF][START_REF] Prieto Barja | Haplotype selection as an adaptive mechanism in the protozoan pathogen Leishmania donovani[END_REF]. Previous studies allowed us to link these forms of genome instability to fitness gain in vitro as judged by the highly reproducible karyotypic changes observed during culture adaptation in independent clinical and animal-derived L. donovani isolates [START_REF] Dumetz | Modulation of Aneuploidy in Leishmania donovani during Adaptation to Different In Vitro and In Vivo Environments and Its Impact on Gene Expression[END_REF][START_REF] Prieto Barja | Haplotype selection as an adaptive mechanism in the protozoan pathogen Leishmania donovani[END_REF][START_REF] Bussotti | Leishmania Genome Dynamics during Environmental Adaptation Reveal Strain-Specific Differences in Gene Copy Number Variation, Karyotype Instability, and Telomeric Amplification[END_REF]. Positive selection of chromosome amplification is further sustained by the independent evolutionary experiment conducted in this study, which once more revealed amplification of chromosomes 5 and 26 as key drivers of in vitro fitness gain. Such karyotypic changes are not exclusive to culture adaptation but have been documented in L. donovani tissue amastigotes [START_REF] Prieto Barja | Haplotype selection as an adaptive mechanism in the protozoan pathogen Leishmania donovani[END_REF],

and in drug resistant clinical Leishmania isolates [START_REF] Laffitte | Plasticity of the Leishmania genome leading to gene copy number variations and drug resistance[END_REF]. Similar to stress-adaptation in fungi [START_REF] Tsai | A Double-Edged Sword: Aneuploidy is a Prevalent Strategy in Fungal Adaptation[END_REF],

karyotypic changes thus may provide the genetic diversity required for Leishmania to evolve beneficial phenotypes in response to environmental change. However, such structural mutations simultaneously affect the expression level of hundreds of genes, raising questions on the nature of the coding sequences that drive karyotypic selection during parasite adaptation, and on the mechanisms that suppress deleterious gene dosage effects while preserving beneficial ones. Applying an integrative systems approach on promastigote parasites during culture adaptation (early passage, EP and late passage, LP) allowed us to address these important open questions.

Comparative genomic, transcriptomic and proteomic analyses of EP and LP promastigote parasites revealed a gene dosage-dependent increase in mRNA and protein abundance for genes implicated in RNA turnover, including RNA-binding proteins known to regulate mRNA stability [START_REF] Haile | Developmental regulation of gene expression in trypanosomatid parasitic protozoa[END_REF],

pumillo domain proteins known to regulate ncRNA abundance [START_REF] Goldstrohm | Post-transcriptional Regulatory Functions of Mammalian Pumilio Proteins[END_REF], and a series of proteins that were associated with trypanosomatid mRNA-binding and post-translational regulation in recent, genome-wide functional screens [START_REF] Erben | A genome-wide tethering screen reveals novel potential posttranscriptional regulators in Trypanosoma brucei[END_REF][START_REF] Lueong | Gene expression regulatory networks in Trypanosoma brucei: insights into the role of the mRNA-binding proteome[END_REF][START_REF] De Pablos | The mRNA-bound Proteome of Leishmania mexicana: Novel Genetic Insight into an Ancient Parasite[END_REF]. Positive selection of chromosome amplifications during L.

donovani culture adaptation thus is likely driven by genes that establish an adaptive, posttranscriptional interface. This interface may regulate differential mRNA stability during fitness gain, which can compensate for deleterious gene dosage effects by selective mRNA degradation, while at the same time boosting the expression of beneficial genes conferring stability to selected mRNAs.

Assessment of gene-dosage-independent expression changes indeed correlated both mRNA degradation and stabilization to the observed fitness phenotype. A number of transcripts implicated in flagellar biogenesis showed reduced stability during fitness gain, which correlated with reduced mobility (data not shown). This coordinated process likely involves shared cis-regulatory sequence elements in the transcripts' 3'UTR that are recognized by the RNA-binding proteins [START_REF] Clayton | Gene expression in Kinetoplastids[END_REF]. In culture, fitness (defined as reproductive success of a given population) is largely synonymous to proliferation, which depends on the number of ribosomes and the cell's translational potential [START_REF] Thomas | An encore for ribosome biogenesis in the control of cell proliferation[END_REF]. However, fitness-adapted mRNA turnover not only changed the quantity but also the quality of the Leishmania ribosomes as judged by the differential expression of various 40S and 60S ribosomal protein isoforms in EP and LP parasites (see Supplementary Tables 2 and5). Such dynamic regulation of ribosomal biogenesis suggests the formation of specialized ribosomes, which not only may increase translation efficiency in these fast-growing LP parasites, but also control translation of unwanted mRNAs, thus providing an additional filter (next to differential RNA stability) to eliminate toxic gene dosage effects. Changes in expression and modification of different rRNA genes, ribosomal proteins, and translation factors indeed can control preferential translation of different subsets of mRNAs in other organisms [START_REF] Guo | Specialized ribosomes and the control of translation[END_REF], including stage-specific translation in Dictyostelium discoideum [START_REF] Ramagopal | Regulation of synthesis of cell-specific ribosomal proteins during differentiation of Dictyostelium discoideum[END_REF][START_REF] Ramagopal | Induction of cell-specific ribosomal proteins in aggregation-competent nonmorphogenetic Dictyostelium discoideum[END_REF], zebrafish development [START_REF] Locati | Expression of distinct maternal and somatic 5.8S, 18S, and 28S rRNA types during zebrafish development[END_REF], or cancer [START_REF] Belin | Dysregulation of ribosome biogenesis and translational capacity is associated with tumor progression of human breast cancer cells[END_REF]. Conducting a dedicated analysis of the L. donovani non-coding (nc) RNome, we have provided further support to the existence of such fitness-adapted ribosomes in Leishmania. First, we observed post-transcriptional upregulation of a large number of snoRNAs and five (out of a total of 9) pseudouridylate synthases in LP compared to EP promastigotes (see Supplementary Table 4). Second, these snoRNA expression changes correlated to changes in the pseudouridinylation (Ψ) profile of the rRNA peptidyl transferase center (PTC) that catalyzes peptide bond formation and peptide release [START_REF] Polacek | The ribosomal peptidyl transferase center: structure, function, evolution, inhibition[END_REF]. Similar Ψ hyper-modification of rRNA was previously observed in bloodstream form trypanosomes and likely contributes to stage-specific adaptation [START_REF] Chikne | A pseudouridylation switch in rRNA is implicated in ribosome function during the life cycle of Trypanosoma brucei[END_REF].

Given the high coding density of chr 5 and chr 26 for snoRNAs and the functional enrichment of these chromosomes for GO term 'rRNA processing' (see Figure S4C), it is interesting to speculate that their convergent amplification in all our evolution experiments may be driven by their ncRNA content and their requirement for fitness-adapted translation. Indeed, defects in rRNA pseudouridylation affect ribosomal ligand binding and translational fidelity in eukaryotic cells [START_REF] Jack | rRNA pseudouridylation defects affect ribosomal ligand binding and translational fidelity from yeast to human cells[END_REF], and changes in PTC modification were shown to affect both the ribosome structure and activity in yeast [START_REF] King | Ribosome structure and activity are altered in cells lacking snoRNPs that form pseudouridines in the peptidyl transferase center[END_REF]. The combination of (i) different rRNA isoforms, (ii) hundreds of snoRNAs and differentially modified rRNA sites, (iii) diverse 40S and 60S ribosomal proteins, and (iv) the formation of different translation complexes [START_REF] Meleppattu | The eIF3 complex of Leishmania-subunit composition and mode of recruitment to different cap-binding complexes[END_REF][START_REF] Yoffe | Evolutionary changes in the Leishmania eIF4F complex involve variations in the eIF4E-eIF4G interactions[END_REF][START_REF] Yoffe | Binding specificities and potential roles of isoforms of eukaryotic initiation factor 4E in Leishmania[END_REF] defines a vast ribosomal landscape in Leishmania. Translational control via fitnessadapted ribosomes likely fine-tunes expression and provides proteomic and phenotypic robustness to adapting parasite populations, which thus can maintain genetic diversity and evolvability despite constant natural selection [START_REF] Prieto Barja | Haplotype selection as an adaptive mechanism in the protozoan pathogen Leishmania donovani[END_REF].

In conclusion, our data uncover Leishmania evolutionary adaptation as an emergent property of a highly complex process that integrates variations in gene dosage with correlating changes in transcript abundance for genes implicated in post-transcriptional regulation and the establishment of specialized ribosomes capable to compensate for toxic gene dosage effects (see Figure 5G). Even though our results are largely correlative in nature, our model is supported by the convergence of the genomic, transcriptomic and proteomic signals we observed between independent populations, which thus are the result of natural selection rather than random genetic drift. Our findings challenge the current genome-centric approach to Leishmania epidemiology and identify the Leishmania noncoding RNome as well as regulatory circuits at transcriptional and translational levels as potential novel sources for biomarker discovery. Finally, our model may be of relevance to other pathogenic systems that gain fitness through genome instability, including fungal infection and cancer.

  animal facilities (A3) accredited by the French Ministry of Agriculture for performing experiments on live rodents (agreement A75-15-01). 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 (dha190013 and 180091) and by the Ministère de l'Enseignement Supérieur, de la Recherche et de l'Innovation (project n°#19683).

  next assessed this phenomenon at genome-wide level by normalizing the RNA-seq read counts to the corresponding DNA-seq reads. Direct comparison of the normalized transcript output in EP.1 versus LP.1 revealed a global, gene dosage-independent increase in transcript abundance in LP.1 log, as would be expected for the increased metabolism in these fast-growing cells (Figure 3C, left panel). In the absence of transcriptional regulation in Leishmania, such a global increase may indicate important changes in RNA turnover, i.e. increased mRNA stability. No difference was observed for metabolically less active, cell-cycle arrested EP.1 and LP.1 stat (Figure 3C, right panel).

  Loss of flagellar function during L. donovani fitness gain has been observed in other independent evolutionary experiments (Bussotti et al., in preparation) and thus represents a reproducible, convergent phenomenon that may liberate ATP for energetically demanding processes that are under positive selection during culture adaptation. Indeed, a large number of transcripts implicated in highly energy-demanding ribosomal biogenesis and translation were stabilized in LP parasites. The differential regulation of mRNA abundance observed in our experimental evolution system thus establishes a first link of post-transcriptional regulation to Leishmania fitness gain.

Hubert

  Curien Maimonide 2018 (GFS and SM), the seeding grant from the Institut Pasteur International Department to the LeiSHield Consortium and the EU H2020 project LeiSHield-MATI -REP-778298-1 (GB), 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 (HV, RL, CP), and the ERD Funds, project CePaViP (CZ.02.1.01/0.0/0.0/ 16_019/0000759) (BV, JS, PV). We thank the CEA-CNRGH for its contribution to the sequencing costs and all the CEA-CNRGH staff who contributed to sample preparation and sequencing for their excellent technical assistance.

Figure 1 :

 1 Figure 1: Phenotypic analysis of EP.1 and LP.1 parasites reveals fitness tradeoff during culture adaptation.

Figure 2 :

 2 Figure 2: RNA-seq analyses of EP.1 and LP.1 promastigotes reveal stage-specific changes in RNA abundance and RNA signatures of fitness gain and fitness cost.

Figure 3 :

 3 Figure 3: RNA abundance during fitness gain is regulated by gene dosage and post-transcriptional mechanisms.

Figure 4 :

 4 Figure 4: Quantitative analysis of the fitness-adapted proteome.

Figure 5 :

 5 Figure 5: The fitness tradeoff in LP promastigotes correlates with snoRNA expression changes and increased rRNA pseudouridinylation levels.

Figure S5 :

 S5 Figure S5: GO analysis of gene dosage-dependent and -independent changes in RNA abundance.

Figure S6 :

 S6 Figure S6: Quantitative proteomics analysis. (A) Box plots representing the median ratio and the upper and lower quartiles of the LFQ intensity values for all LP biological replicates (LP.2, LP.3, LP.4 and LP.5) compared to the median of all EP replicates (EP.2, EP.3, EP.4 and EP.5). (B) Cluster analysis of all EP and LP samples (Ward method). (C) Ratio plot representing the mean LFQ intensity value between EP and LP for each individual, quantified protein. The experimental and the expected regression lines are shown in red and blue respectively. (D, E) Cluster efficiency for the GO category 'biological process' for proteins whose abundance correlates with increased (D) or decreased normalized DNA read counts (E) in LP log parasites. (F) Cluster efficiency for genes showing gene dosage-independent changes in protein abundance. (G) Table listing selected genes associated with the GO term 'ciliary cell motility' from the GO enrichment analysis presented in Figure 4E. Their respective fold change values computed from Protein LFQ intensities (grey bars) and DNA normalized read counts (black bars) for LP versus EP log parasites are represented.

Figure S7 :

 S7 Figure S7: Enrichment of small RNAs obtained from post-ribosomal supernatants (PRS) of EP promastigotes. 2x10 9 cells were disrupted by nitrogen cavitation under low salt concentration (150 mM KCL) in the presence of high MgCl 2 (10 mM) followed by ribosome extraction using high KCL (300 mM). The ribosomes were removed by centrifugation at 35,000 rpm for 2h. 2µg of RNA from total lysate (Total RNA) and PRS sample were separated on a 10% polyacrylamide gel and stained with ethidium bromide.

  Figure S7

  1 and LP.1 log, stat, and EP.1 meta parasites. The low yield in LP.1 meta parasites precluded their analysis by RNA-seq.

	Principal component and hierarchical clustering analyses demonstrated that transcript profiles of
	EP.1 and LP.1 parasites grouped according to stage, indicating that stage-specific expression changes
	in log, stat and meta parasites dominate over those associated with the EP.1/LP.1 promastigote
	fitness tradeoff (Supplementary Table 2a-f, Figure 2A, Figure S3A). Significant stage-specific changes
	were observed in EP.1 and LP.1 parasites during the log-stat transition for respectively 54.2% and
	49.3% of the transcripts and ca. 35% of the promastigote transcriptome was modulated between
	EP.1 stat and EP.1 meta (Figure 2B and S3B, Supplementary Table 2a-f).

We next assessed changes in transcript abundance observed at logarithmic growth phase in LP.1 compared to EP.1 promastigotes to gain first insight into pathways associated with in vitro fitness gain (i.e. accelerated growth). We identified 345 transcripts with increased abundance in LP.1 log (Figure

Acknowledgements

This work was supported by the Fondation de la Recherche Médicale contract FDT201805005619 (LP), 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, LP, PP, GB), the Campus France Franco-Israeli Programme

Supplementary

*Due to the volume of the supplementary Tables, their submission to bioRxiv as pdf files was deemed unfeasible. All supplementary Tables can be requested from the corresponding authors.