Pablo Prieto Barja

Pascale Pescher

Giovanni Bussotti

Franck Dumetz

Hideo Imamura

Darek Kedra

Malgorzata Domagalska

Victor Chaumeau

Heinz Himmelbauer

Michel Pagès

Michel Pages

Yvon Sterkers

Jean-Claude Dujardin

Cedric Notredame
email: cedric.notredame@crg.eu

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

Haplotype selection as an adaptive mechanism in the protozoan pathogen Leishmania donovani

Keywords: Leishmania donovani, aneuploidy, genome instability, dosage compensation, fitness, haplotypes

Introduction

Rapid pathogen adaptation to novel environments is a major threat to human health. In parasites evolving under constant host selection, fitness gains often cause increased pathogenicity [START_REF] Pallen | Bacterial pathogenomics[END_REF] . Aneuploidy recently emerged as a key driver for evolutionary adaptation of fungal and protist pathogens [START_REF] Torres | Effects of aneuploidy on cellular physiology and cell division in haploid yeast[END_REF][START_REF] Pavelka | Aneuploidy confers quantitative proteome changes and phenotypic variation in budding yeast[END_REF][START_REF] Pavelka | Dr Jekyll and Mr Hyde: role of aneuploidy in cellular adaptation and cancer[END_REF] , in which chromosome copy number variations induce phenotypic changes through modulation of transcript and protein abundance [START_REF] Selmecki | Acquisition of aneuploidy provides increased fitness during the evolution of antifungal drug resistance[END_REF] . In addition to its effect on gene count, aneuploidy impacts cellular phenotypes through the selection of beneficial alleles. While this phenomenon has received considerable attention in human genetics [START_REF] Sabeti | Detecting recent positive selection in the human genome from haplotype structure[END_REF] and cancer genome analysis [START_REF] Staaf | Landscape of somatic allelic imbalances and copy number alterations in human lung carcinoma[END_REF][START_REF] Scheinfeldt | Recent human adaptation: genomic approaches, interpretation and insights[END_REF][START_REF] Gao | Punctuated copy number evolution and clonal stasis in triple-negative breast cancer[END_REF] , the role of aneuploidy on allelic selection in genome evolution and adaptation in pathogenic eukaryotes remains poorly understood in most organisms, in spite of recent attention [START_REF] Sterkers | Parasexuality and mosaic aneuploidy in Leishmania: alternative genetics[END_REF][START_REF] Imamura | Evolutionary genomics of epidemic visceral leishmaniasis in the Indian subcontinent[END_REF] .

We addressed this question in the protozoan parasite Leishmania donovani, an important human pathogen that causes fatal visceral leishmaniasis [START_REF] Alvar | Leishmaniasis worldwide and global estimates of its incidence[END_REF] . During its life cycle Leishmania undergoes a major developmental transition from insect-stage promastigotes to mammalian-stage amastigotes, which adapts these parasites for extra-and intracellular survival, respectively [START_REF] Zilberstein | The role of pH and temperature in the development of Leishmania parasites[END_REF] . In addition to environmentally-induced stage differentiation, Leishmania can adapt to a variety of unpredictable fluctuations inside its human host, notably pharmacological interventions. Such environment-genotype interactions likely select parasites for higher fitness and have important consequences on disease outcome as demonstrated by the emergence of drug resistant clinical isolates [START_REF] Leprohon | Drug resistance analysis by next generation sequencing in Leishmania[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] . However, the underlying genetic mechanisms that drive short-term Leishmania evolution are still largely unknown. In the absence of classical transcriptional gene regulation, these early-branching eukaryotes often control protein abundance via gene and chromosome amplification [START_REF] Sterkers | Parasexuality and mosaic aneuploidy in Leishmania: alternative genetics[END_REF][START_REF] Dujardin | Mosaic aneuploidy in Leishmania: the perspective of whole genome sequencing[END_REF] . This process is fueled by frequent asymmetric chromosome allotments, which cause constitutive mosaicism [START_REF] Rogers | Chromosome and gene copy number variation allow major structural change between species and strains of Leishmania[END_REF][START_REF] Sterkers | FISH analysis reveals aneuploidy and continual generation of chromosomal mosaicism in Leishmania major[END_REF] . This model has recently received new support [START_REF] Imamura | Evolutionary genomics of epidemic visceral leishmaniasis in the Indian subcontinent[END_REF] with a largescale deep sequencing study of L. donovani field isolates showing a correlation between polysomy variations and changes in heterozygosity. Leishmania therefore represents an ideal non-conventional system to investigate how karyotype variations and allelic selection drive parasite fitness gains and how transient and stable aneuploidies impact the long-term evolutionary trajectory.

Results

Karyotype and haplotype selection in field isolates

We first explored the impact of aneuploidy on allelic selection in Leishmania using the recently published sequencing data of 204 L. donovani field isolates [START_REF] Imamura | Evolutionary genomics of epidemic visceral leishmaniasis in the Indian subcontinent[END_REF] . These populations originate from independent evolutionary radiations that occurred after a population bottleneck caused by DDT vector control campaigns in the 1960s. This makes them a convenient benchmark to study genetic diversity dynamics. We conducted read-depth analyses on publicly available sequencing data sets [START_REF] Imamura | Evolutionary genomics of epidemic visceral leishmaniasis in the Indian subcontinent[END_REF] to explore ploidy variations across culture-adapted isolates (Fig. 1a). Even though a large number of intermediate read-depth values suggest frequent mosaicism (Fig. 1a, Supplementary Table 1), the relatively high frequency of discrete polysomy levels other than disomic levels (notably trisomy) indicates that a large fraction of isolates may be homogenous for single or co-occurring aneuploidies (e.g chr 5, 9, 23 and 26, Fig. 1b). While all chromosomes are predominantly disomic except for chr 31, the overall karyotypic diversity is very high with the most frequent karyotype representing only about 10% of the 204 isolates.

The frequent occurrence of karyotypic similarities among isolates suggests that the corresponding aneuploidies may either have a common, monoclonal origin (i.e. one unique founding duplication) or result from frequent convergent processes (i.e. selection of several independent duplications of the same homologous chromosome). Relatively weak correlations between aneuploidy and SNP variation [START_REF] Imamura | Evolutionary genomics of epidemic visceral leishmaniasis in the Indian subcontinent[END_REF] does not discriminate between these two possibilities. While clonal expansion of a single founder may underlie the most common aneuploidies (like chr 31, 5 or 26), their occurrence in karyotypically very diverse isolates is also compatible with a polyclonal origin.

Under a polyclonal scenario, each culture-adapted isolate arises from the expansion of mixed subpopulations sharing some independently generated polysomies. This model is in good agreement with mosaicism, previously observed in cultured Leishmania major [START_REF] Sterkers | FISH analysis reveals aneuploidy and continual generation of chromosomal mosaicism in Leishmania major[END_REF] . We analyzed allelic frequency distributions to evaluate the likelihood of polyclonal origin (Supplementary Fig. 1). In this analysis, well-defined peaks in the merged profiles of most chromosome allelic frequency distributions indicate that most heterozygote alleles are balanced similarly across isolates. For disomic profiles, we typically observed unimodal allelic frequency distributions centered at 50% (Supplementary Fig. 1a -left column). These observations are consistent with diploid asexual populations in which most alleles diffuse by random drift.

Chromosomes with a read-depth level compatible with trisomy present us with a more complex situation. In these chromosomes, the merged allelic frequency profiles (Supplementary Fig. 1a -right column) revealed a continuum between the two extreme cases of pure unimodal (e.g. chr 1, Fig. 1c, Supplementary Fig. 1d bottom scenario) and pure bimodal (e.g. chr 5, Fig. 1c, Supplementary Fig. 1d other scenarios) allelic frequency distributions. In the absence of mosaicism, trisomies usually cause bimodal allelic frequency distributions. Unimodal distributions like chr 1 therefore reflect population mosaicism as a consequence of a duplication scenario with both homologous chromosomes equally likely to be duplicated across the entire population of cells constituting an isolate, thus maintaining the original allelic frequency equilibrium of the disomic state (Supplementary Fig 1d). The coexistence of haplotypically different trisomies within the same isolate clearly indicates mosaicism, which strongly suggests a polyclonal origin. In contrast, bimodal allelic frequency distributions centered on 33 and 66%, like that observed for chr 5, reflect the collective doubling in relative frequency of half of the alleles and are characteristic of homogenous duplication of the same homologous chromosome within each isolate population (thus defining an identical haplotype). While the integration of various mosaicism combinations by sequencing analysis could cause similar read-depths, it is unlikely that mosaicism would also cause bimodal allelic frequency distributions centered on 33 and 66%.

Even though they occur in a large fraction of the isolates, chr 5 trisomies are most likely the result of a convergent process as suggested by their polyphyletic dispersion on the phylogenetic tree of the 204 filed isolates [START_REF] Imamura | Evolutionary genomics of epidemic visceral leishmaniasis in the Indian subcontinent[END_REF] . Given such a convergent process, we asked if these trisomies were selected as karyotypes or as specific haplotype combinations.

Haplotype selection should cause increased homogeneity across isolates. We therefore measured allelic dispersion across trisomic isolates. We estimated the relative loss of heterozygosity (LOH) between disomic and trisomic isolates for each chromosome. This analysis shows that chr 5 had one of the lowest allelic diversity (highest average LOH) and suggests frequent selection of the same homologous chromosome among the trisomic isolates.

On these same subsets of trisomic and disomic chromosomes, we measured the segregation variation (SV). SV is defined as the average ratio between the frequency of the dominant allele in each trisomic isolate and its corresponding average frequency in disomic isolates. This index reflects the changes in segregation patterns between disomic and trisomic states and is highest for chromosomes that show a unimodal distribution at a disomic state and a bi-modal distribution at a trisomic state. As expected, we found chr 5 to have one of the highest SV ratio. We also found a strong correlation between LOH and SV across all chromosomes with the scatter plot displaying a clear continuum between chr 1 and chr 29 (Fig. 1d). This observation is in agreement with our hypothesis that high levels of haplotype selection cause increased allelic homogeneity across isolates. These results confirm that bimodal allelic frequency profiles are more likely induced by a selection process rather than a monoclonal founding event. Monoclonal origin without haplotype selection would cause a high SV and low LOH. Interestingly, distinctly bimodal chr 5 trisomies coexist with unimodal trisomies of other chromosomes in a subset of the isolates (e.g. chr 12, Supplementary Fig. 2). This finding indicates that various levels of haplotypic selection can affect different chromosomes within the same parasite isolate. Altogether, the complex karyotypic and haplotypic variations across field isolates clearly favors a polyclonal mosaic origin as the most common scenario with diversity drawn from an equally diverse original population, as postulated by Sterkers et al. [START_REF] Sterkers | FISH analysis reveals aneuploidy and continual generation of chromosomal mosaicism in Leishmania major[END_REF] .

In vivo karyotype fluctuations

Polyclonal origin implies a pre-existing population-wide karyotypic mosaicism within each individual isolate. Maintaining this degree of variability would require frequent polysomic fluctuations. We measured this effect by monitoring chromosomal copy number when repassaging an aneuploid field isolate (strain LdBPK282) through hamsters. The rapid shift to a predominant disomic karyotype in the field isolate analysed here (and an independent isolate shown by Dumetz et al. [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] ) demonstrates highly dynamic aneuploidy fluctuations (Fig. 2a). It is unclear, however, if the observed aneuploidies during culture adaptation arise from expansion of existing subpopulations or constitute a reversible de novo phenomenon. We therefore examined this question applying single cell DNA-FISH analysis on liver and spleen amastigotes. We analyzed the experimental Sudanese L. donovani strain LD1S with available DNA-FISH probes to monitor aneuploidies for chr 5, 17, 22 and 27. Hamsters were infected with amastigotes isolated from infected hamster spleens that showed a disomic read depth for all chromosomes except tetrasomic chr 31 (see Supplementary Fig. 3a). For all chromosomes investigated, our analysis confirmed frequent in situ mosaic aneuploidies in liver and spleen (Fig. 2b). This pre-existing diversity can facilitate population-wide aneuploidies in response to environmental changes. This possibility is in line with the mosaic origin of aneuploidies observed in culture [START_REF] Sterkers | FISH analysis reveals aneuploidy and continual generation of chromosomal mosaicism in Leishmania major[END_REF] .

In vitro karyotype and haplotype selection

While field isolate and in vivo analyses suggest that parasite adaptation relies on a mechanism dedicated to accommodate frequent aneuploidies, these observations only provide static snapshots of apparently highly dynamic fluctuations. To further elucidate this process, we used again the experimental strain LD1S. In contrast to L. donovani clinical isolates from the Indian sub-continent, this strain contains a larger number of heterozygous sites (23014 compared to a median of 1143 per field isolate), which makes it an ideal model for longitudinal monitoring of chromosome copy number variation and haplotype selection.

We used HTSeq to follow hamster-derived LD1S amastigotes during adaptation to in vitro culture. Read depth analyses indicated that parasites rapidly establish stable trisomies for chr 5, 9, 23, and 26 between in vitro passages p2 (ca. 20 generations) and p10 (ca. 60 generations) (Fig. 3a and Supplementary Fig. 3a). These karyotypic trajectories matched the most common variations observed in the field isolates (Fig. 1a andb, Supplementary Fig. 1) and they were highly reproducible across two independent experiments [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] . Not all aneuploidies are stable and homogenous, however. For example, chr 20 underwent a transient trisomy between passages p2 and p20 (ca. 190 generations), while probable mosaic aneuploidies (Fig. 3a) occurred for chr 14 and 15 and were stably maintained as judged by their intermediate read-depth levels at p10 and p20.

Population mosaicism in Leishmania major was previously examined on cultured promastigotes using a combination of sub-cloning and DNA-FISH analysis [START_REF] Sterkers | FISH analysis reveals aneuploidy and continual generation of chromosomal mosaicism in Leishmania major[END_REF] . We used systematic comparison between p20 and eight individual sub-clones to model the original population complexity (Fig. 3b and Supplementary Fig. 3b). Haplotype and karyotype comparisons suggest that the eight clones may have arisen from at least two independent founding individuals (Supplementary Fig. 3c). Indeed, aneuploidies provide a powerful insight into haplotype phasing, due to systematic frequency shifts associated with duplicated homolog chromosome [START_REF] Legrand | Haplotype mapping of a diploid non-meiotic organism using existing and induced aneuploidies[END_REF] . In this case, direct haplotype comparisons (Fig. 4a, Supplementary Fig. 4) clearly show how the allelic nature of chr 5 and chr 9 trisomies sets clones 1 and 8 apart from the rest. For instance, the existence of two groups with differences of major alleles on chr 5 (1 and 8 vs the rest) indicates that these trisomies were established through duplications of different chr 5 homologs (i.e. same chr 5 homologue duplicated in clones 1 and 8 that is different from the other clones). The most parsimonious reconstruction simply requires two trisomies of independent origin and is consistent with a high aneuploidy turnover [START_REF] Sterkers | FISH analysis reveals aneuploidy and continual generation of chromosomal mosaicism in Leishmania major[END_REF] . This explanation also applies to chromosome 9 and further supports the common origin of clones 1 and 8. Sequencing 8 subclones has limited statistical power, and our speculation of two original founders does rely on a parsimonious hypothesis. Yet, the existence of clearly distinct subpopulations featuring similar aneuploidies further supports the mosaic origin scenario, under which individual strains represent complex mixtures of stable, karyotypically distinct subpopulations.

Between subpopulations, similar karyotypes and haplotypes likely reflect a form of convergence. At the karyotypic level, the most obvious indications of convergence occur on chr 5, 9 and 26, whose trisomies are shared across all clones despite their independent origins. Selection could also occur at the haplotypic level, with the clearest indication provided by chr 26 whose haplotype map shows a consistent selection of the same allelic combination across clones of different origin (Fig 4a). Unfortunately, the simultaneous haplotype/karyotype selection taking place on this chromosome makes it impossible to discriminate between purifying (lethality of one of the two possible trisomies) or positive (higher fitness of one trisomy) selection pressure being imposed on its haplotype.

Chromosome 20 instead provides clear evidence for positive haplotype selection. This chromosome is disomic and heterozygous at p2, becomes trisomic at p10, but reverts back to disomy at p20 (Figs. 3a and4b, Supplementary Fig. 5). After reversion, the allelic frequency profiles unexpectedly show a near-perfect homozygote disomy in 6 out of 8 clones (i.e. relatively flat allele frequency profile combined with a 60% decrease in the number of heterozygous sites, Fig. 4a and Supplementary Fig. 6). In chr 20, transient trisomy therefore provides an intermediate step to establish a homozygous disomy, whose selection contributes to positive fitness as previously seen on yeast [START_REF] Yona | Chromosomal duplication is a transient evolutionary solution to stress[END_REF] .

Correlation of aneuploidies with phenotypic variations

Since Leishmania primarily lacks classical gene regulation mechanisms, it has long been speculated that polysomy should result in proportional transcriptional fluctuations [START_REF] Downing | Whole genome sequencing of multiple Leishmania donovani clinical isolates provides insights into population structure and mechanisms of drug resistance[END_REF] and that aneuploidy mediated gene expression changes could be a trait under selection. We used variation across each individual isolate shows that polysomy-prone chromosomes exhibit a significantly higher level of heterozygous sites than their more stable counterparts. This trend is especially strong for chr 31, which has the most stable, frequent and highest order aneuploidy. This finding suggests that L. donovani might utilize aneuploidy to accumulate mutations and increase its diversity. Further, these results confirm that even though transient in vivo aneuploidies are difficult to detect and quantify, they are frequent enough to shape the parasite genome and its genetic diversity.

Discussion

Utilizing recently published genomes of 204 field isolates [START_REF] Imamura | Evolutionary genomics of epidemic visceral leishmaniasis in the Indian subcontinent[END_REF] and conducting evolutionary experiments, we demonstrated that karyotypic variation modulates transcript abundance and thus likely generates considerable phenotypic variability with possible fitness gains associated with tissue-specific haplotype selection in situ. The genomic landscape defined by an exhaustive combination of all possible trisomies, disomies and monosomies along with their haplotype variations is substantial and at least one order of magnitude larger than a similar estimate recently made on the basis of aneuploidy variations [START_REF] Sterkers | Novel insights into genome plasticity in Eukaryotes: mosaic aneuploidy in Leishmania[END_REF] . Aneuploidy turnover therefore provides an efficient alternative to sex-based haplotype selection. This mechanism may be especially useful within the mammalian host, where the absence of sexual reproduction creates an even heavier selective pressure on the parasite's survival prospects.

Rapid aneuploidy turnover combined with haplotype selection permits fast adaptation but comes with a significant genetic cost. Since alleles with lowered frequencies are at a higher risk of disappearing from the genetic pool [START_REF] Lande | Natural Selection and Random Genetic Drift in Phenotypic Evolution[END_REF] , haplotype selection may induce loss of heterozygosity. The transient trisomy of chr 20 and its in vivo, tissue-specific haplotypic diversity illustrates the dilemma for the parasite between over-adaptation to any given environment that may involve irreversible genetic tradeoffs and maintenance of sufficient genetic diversity for future adaptations. In this chromosome, the loss of heterozygosity could contribute to adaptation across conditions and tissues, however this same process may also yield an evolutionary dead end. The evolutionary success of Leishmania testifies of its ability to establish a delicate balance between these conflicting requirements for short-and longterm adaptation. Long-term survival for the parasite involves two complementary mechanisms to generate and maintain genetic diversity. The first one relies on the mosaic origin of selected populations. We demonstrated that culture adapted isolates do not have a single founding parent but result from the simultaneous selection of several individuals. This process maintains a high level of karyotypic and haplotypic diversity allowing polysomies to occur frequently and independently in genetically distinct individuals. The second source of diversity is a direct consequence of the relaxed selection that occurs after gene or chromosome duplication [START_REF] Lynch | The evolutionary fate and consequences of duplicate genes[END_REF] . Our findings agree with these models and further suggest that aneuploidy-prone chromosomes have significantly higher mutation rates than their more stable counterparts. Thus, aneuploidy represents an additional powerful mechanism to increase genetic diversity in L. donovani.

Leishmania, like all known infectious agents, provides us with an excellent example of genetic selection driven by survival. But it does so in a very unusual way: while high mutation rates and genetic material exchange through sexual or analogous mechanisms mediate survival in most pathogens, Leishmania evolved an additional strategy based on genome instability to enhance parasite evolvability [START_REF] Victoir | How to succeed in parasitic life without sex? Asking Leishmania[END_REF][START_REF] Bastien | Leishmania: sex, lies and karyotype[END_REF][START_REF] Dujardin | Clonal propagation and the fast generation of karyotype diversity: An in vitro Leishmania model[END_REF] . We show here that the parasite harnessed this process and efficiently combined it with haplotype selection for rapid and efficient genome adaptation under environmental pressure. This strategy has direct consequences for current protocols in Leishmania drug and biomarker discovery. First, parasite biology and select for resistance phenotypes [START_REF] Gazanion | Cos-Seq for high-throughput identification of drug target and resistance mechanisms in the protozoan parasite Leishmania[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] Leprohon | Drug resistance analysis by next generation sequencing in Leishmania[END_REF] . New strategies for antileishmanial drug discovery should avoid direct parasite selection and rather target the parasites' dependence on host cell metabolism. Second, the massive genomic changes we observed during Leishmania culture adaptation challenges current protocols for biomarker discovery, which all rely on in vitro expansion of clinical isolates. In the light of our findings, we should modify these protocols to identify additional biomarkers hidden/diluted by culturebased methods especially when polysomy itself is considered a biomarker. Alternatively, we could establish culture-independent approaches that propagate field isolates in experimental animal models or apply direct tissue sequencing. It is clear that genome instability needs to be considered when investigating medically relevant aspects of the parasite's phenotype, such as tissue tropism, drug susceptibility and pathogenicity. Our findings establish a foundation for the future discovery of haplotypes with diagnostic and prognostic value.

Material and Methods

L. donovani strains, culture conditions and cell cloning. Culture-adapted L. donovani field isolates maintained for more than 20 in vitro passages from the Indian Sub-Continent were used in our analyses to track evolutionary diversity maintained across populations [START_REF] Imamura | Evolutionary genomics of epidemic visceral leishmaniasis in the Indian subcontinent[END_REF] .

Isolates were not sub-cloned and populations were used for all studies. Infectious L. donovani strain 1S2D (MHOM/SD/62/1S-CL2D, referred to here as LD1S) was obtained from Henry Murray (Weill Cornell Medical College, New York, USA) and maintained in hamsters by serial passaging (see below). All LD1S studies were performed with populations except when specified otherwise. For promastigote differentiation and culture, 1 x 10 7 LD1S amastigotes purified from hamster spleens were cultured at 26 °C in M199 media supplemented with 10 % FCS, 25 mM HEPES pH 6.9, 4 mM NaHCO 3 1 mM glutamine, 1 x RPMI 1640 vitamin mix, 0.2 μM folic acid, 100 μM adenine, 7.6 mM hemin, 8 µM biopterin, 50 U/ml of penicillin, and 50 μg/ml of streptomycin. Promastigotes were then maintained in culture by serial dilution in fresh medium once they reached stationary phase. At passages 2, 10 and 20 corresponding to approximately 20, 60 and 190 generations respectively, parasites in exponential growth phase were collected and adjusted to 2 x 10 8 parasites per tube for DNA extraction (see below). After 20 in vitro passages, serial dilutions of promastigotes were plated on M199 Agar plates and 8 clones were selected and amplified in promastigote medium.

Hamster infection and isolation of infectious amastigotes. We used female

RjHan:AURA golden syrian hamsters between 4 to 6 weeks of age (Janvier Labs, France) and a weight between 60 to 90g. Housing of the animals and all experiments were conducted in agreement with the project number 2013-0092 approved by the Institut Pasteur Ethics committee in accordance to the European legislation/guidelines EU 2010/63. We used a total of 13 hamsters for DNA-FISH analysis (1), DNA/RNA extraction from purified amastigotes (2), hamster infection (6), monitor tissue parasite burden (3), and comparative HTseq analysis of liver and spleen amastigotes (1). Given the nature of our animal experiments to compare infected versus non-infected animals and use infected animals as a source of biological material, our studies were not blinded or randomized. No animals were excluded from analysis. Anesthetized hamsters were inoculated by intra-cardiac injection with 5 x 10 7 amastigotes obtained from infected hamster spleens. Hamster weight was monitored and animals were euthanized with CO 2 after four months of infection. Spleens and livers were collected, weighed, and homogenized in PBS as described [START_REF] Pescher | Quantitative proteome profiling informs on phenotypic traits that adapt Leishmania donovani for axenic and intracellular proliferation[END_REF] . For amastigote purification, tissues were homogenized in 25 ml PBS supplemented with 2.5 mg/ml saponine using the gentleMACS homogenizer and gentleMACS M tubes from Miltenyi. Homogenates were cleared by centrifugation at 130g for 5 min at 4°C, the supernatants were collected, 1 ml of saponine (25 mg/ml) was added under gentle agitation, and parasites were harvested 5 min later by centrifugation at 2000g for 10 min at 4°C. After two washing steps with PBS, remaining host cell contaminants were removed by Percoll centrifugation. Briefly, parasite were resuspended in 6 ml of 45 % Percoll, and 3 ml were layered on a cushions of 2 ml of 90 % Percoll in 15 ml falcon tubes. After 35 min of centrifugation at 3500g, 15°C, amastigotes were recovered from the interface of the gradient and washed 3 times in PBS (centrifugation at 2000g, 10min, 4°C). Tissue-derived amastigotes were adjusted to 2 x 10 8 parasites per tube for RNA or DNA extraction or inoculated in culture medium for differentiation into promastigotes and further culture. DNA FISH analysis. DNA probes for chromosomes 5, 17, 22, and 27 were prepared as described in [START_REF] Sterkers | FISH analysis reveals aneuploidy and continual generation of chromosomal mosaicism in Leishmania major[END_REF] . Amastigotes were purified from infected hamster livers and spleens (previously inoculated with amastigotes obtained from infected hamster spleens) as described above, immobilized on glass slides, fixed in 4% paraformaldehyde and 4% acetic acid. Fluorescence in situ hybridization was performed according to Sterkers et al. [START_REF] Sterkers | FISH analysis reveals aneuploidy and continual generation of chromosomal mosaicism in Leishmania major[END_REF] .

Leishmania cells were viewed by phase contrast, and fluorescence was visualized using appropriate filters on a ZeissAxioplan 2 microscope with a 100 x objective. Digital images were captured using a Photometrics CoolSnap CDD camera (Roper Scientific) and processed with MetaView (Universal Imaging). Z-Stack image acquisitions (15 planes of 0.25 mm) were systematically performed for each cell analyzed using a Piezo controller, allowing to view the nucleus in all planes and to count the total number of labeled chromosomes. The ploidy was estimated on 100 labeled cells. We carefully considered the cell cycle state and excluded cells that were replicating their DNA as judged by the presence of large or duplicated kinetoplasts, which corresponded to less than1% and thus a negligible fraction. This low replication rate was confirmed by our own infection data analysis that allowed us to calculate a generation time of amastigotes in vivo of 10.4 days.

Genomic sequencing. All sequencing analysis were performed using Illumina short-read technology. The sequences for BPK282/0 promastigotes and amastigotes shown in Figure 2a were generated using the Nextera XT library preparation kit (Illumina), the KAPA Library Quantification Kits (KAPA Biosystems), and paired-end sequencing on a HTSeq 1500 sequencer (Illumina). The sequences of LD1S liver and spleen amastigotes (Figure 5c, Supplementary Figure 8 and9) were carried out as a service by the Centro Nacional de Análisis Genómico (CNAG, Barcelona) using the following protocol: 0.3-0.6 micrograms of sheared genomic DNA (Covaris E201) was end-repaired, adenylated and ligated to Illumina specific indexed paired-end adapters (adapter:insert molar ratio 40:1). The DNA library was size selected with AMPure XP beads in order to reach an insert size of 220-550bp and to remove unligated adapters. The final libraries were quantified with the Library Quantification Kit (Kapa Biosystems). On the base of the qPCR quantification the molar concentration sufficient for sequencing run was estimated. Each library was sequenced using TruSeq SBS Kit v3-HS (Illumina Inc.), in paired end mode, 2x101bp, in a fraction of one sequencing lane of an HTSeq2000 flowcell v3 (Illumina Inc.) according to standard Illumina operation procedures with a yield between 3,593 -6,326 Gb for each sample. The rest of the LD1S Leishmania sequences were generated at the Centre for Genomic Regulation (CRG, Barcelona) using the following protocol: DNA was isolated using DNeasy blood and tissue kits from Qiagen, nucleic acid concentration was measured with a NanoDrop® spectrometer, and between 2 to 5 µg were used for sequencing and quality control. DNA sequencing was performed using an Illumina HiSeq 2000 platform and TruSeq v3 kits. One µg of nucleic acids was used for DNAseq while the rest of the material was used for quality control. All libraries were sequenced in 125 bp reads. A summary of the read numbers and read-length associated with each genome is provided in Supplementary Table 2.

Read depth Sequencing Analysis. All DNA sequences were mapped against L. donovani BPK282A1 from TriTrypDBv7. Mapping and post processing was carried out using a combination of BWA mem (v0.7.8) [START_REF] Li | Fast and accurate short read alignment with Burrows-Wheeler transform[END_REF] along with Samtools (v1.3) [START_REF] Li | The Sequence Alignment/Map format and SAMtools[END_REF] in order to refine read information and clean mapped sequences. Alignments were further refined using GATK (v2.8) IndelRealigner and MarkDuplicates [START_REF] Mckenna | The Genome Analysis Toolkit: a MapReduce framework for analyzing next-generation DNA sequencing data[END_REF] . Samtools depth command was used to estimate read-depth for every base, so as to determine median read-depth for each chromosome. In order to estimate chromosome polysomy level in each sample, chromosome 34 -which showed a stable disomic level across samples -was used to normalize other chromosomes read depth. This measure was calibrated to determine polysomy levels, under the assumption that the most frequent polysomy levels correspond to disomy. Intermediate polysomy levels were classified as trisomic or tetrasomic events by using cutoffs on the normalized chromosome read depth values. Chromosomes with values ranging between 2.6 and 3.4 were classified as trisomic, and values between 3.6 and 4.4 as tetrasomic. These calibrated values were used to select subsets of chromosomes having comparable aneuploidy levels. (Accession number: PRJEA61817) Aneuploidy co-occurrence. Aneuploidy co-occurrence analysis was carried out by measuring the Pearson correlations between pairs of chromosomes, with each chromosome being represented as vector of normalized depth values, each value being drawn from an isolate. P-Values were measured so as to identify pairs having the highest levels of covariation -as opposed to constant or unrelated polysomy levels.

Variant calling and allele frequency analysis. Pileups were generated across the reference genome for all the samples. Samtools mpileup and its multiallelic caller were then used to detect variants separately in all samples. In order to filter out sites potentially affected by episomal amplifications, regions with a read-depth higher than 1.6 x the median depth of each chromosome were excluded from the analysis. The resulting sites were then used to draw allele frequency statistics, and produce allele frequency profiles. Variants were then filtered using a quality threshold of 15 using a Read Position Bias (RPB) filter of 0.1 and alleles with frequencies comprised between 0.1 and 0.9 were retained. This approach, that was explicitly designed to prioritize sensitivity over specificity, yields about 10 times more sites than a more conservative approach (28071 vs 2418 sites [START_REF] Imamura | Evolutionary genomics of epidemic visceral leishmaniasis in the Indian subcontinent[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] ). The expected rate of false positive associated with this procedure can be shown to be in the order of 3% (-10*log10(1/31.6)). This amount of noise that would be incompatible with biomarker discovery does not compromise the current analysis focused on tracing homologue chromosome fate in aneuploidies.

Merged allele frequency profiles. The allele frequencies associated with variants called as explained in the previous section were then used to build merged allele frequency distribution profiles for each chromosome across isolates. Profiles of homologue chromosomes classified with the same polysomy level -as explained in read depth analysis -were added up to yield a set of merged profiles (Supplementary Figure 1). The stacking of pre-computed profiles presents the advantage of maintaining the individual segregation profiles and can therefore be used to estimate profile homogeneity across samples as indicated by well-defined peaks that reflect conserved allelic frequencies across isolates.

Most of these profiles exhibit the near canonical profiles expected for allele distribution in a diploid, triploid or tetraploid organism and therefore confirm the suitability of the noise to signal ratio of our site calling procedure.

Segregation variations and Loss of Heterozygosity dispersion plots. Segregation variations (SV) between disomic and trisomic chromosomes were estimated by considering

for each heterozygous site the log ratio between the frequency of the most major (major allele) allele in the considered isolate, and the average frequency of the same allele in the disomic population. With s heterozygote site in an isolate t isolate trisomic for chromosome N -T all the isolates trisomic for chromosome N d isolate disomic for chromosome N -D all the isolates disomic for chromosome N -MA(s,t) the Major Allele on site s in the isolate t trisomic for Chromosome N -Frequency (MA(s,t),t) is the frequency of MA(s,t) in isolate t -AvgFrequency (MA(s,t),D) the average frequency of MA(s,t) across all isolates D disomic for Chromosome N This procedure amounts to considering the disomic average as a null hypothesis with respect to the fate of the corresponding alleles in their trisomic states. When doing so the nature of the allele is not considered and two isolates having different major alleles on the same site can nonetheless contribute the same value to the distribution of ratios, as plotted in Supplementary Figure S1b. These ratios were then averaged for each chromosome and used as a measure for Segregation variation. On Supplementary Fig S1, profiles are ordered according to this ratio that tends to reflect the disomic tendency in the trisomics.

Loss of heterozygosity (LOH) was estimated by considering for each trisomic heterozygous site the log odd ratio between the average frequency of the major allele across trisomic isolates and the corresponding frequency across disomic isolates. These ratios were then averaged for each chromosome and used as a loss of heterozygosity measure Haplotype phasing. Haplotype phasing was carried out by considering the effect of aneuploidy on the original frequencies following the protocol originally defined in [START_REF] Legrand | Haplotype mapping of a diploid non-meiotic organism using existing and induced aneuploidies[END_REF] . In this analysis, alleles undergoing similar shifts are considered as being part of the same haplotype block. While heterozygous SNPs are not dense enough in the 204 filed isolates for this approach to be validated, we took advantage of the LD1S hybrid origin to verify that SNPs occurring with the insert distances of the pair-ended sequencing had similar segregation shifts after duplication and could therefore considered to be genetically linked, thus ruling out recombination as a major confounding factor. This finding is in agreement with the original study of very low recombination rates in the 204 filed isolates. Our analysis indicates that 50% of the SNP pairs within insert distance of one another have allelic frequencies less that 10% point apart in the bimodal trisomics (Fig S11 ), as compared with 23% when randomizing pair ends in order to test for pair connectivity and allele linkage.

Transcriptome sequencing. RNA was prepared using RNeasy mini plus blood and tissue kits from Qiagen according to manufacturer instructions. Nucleic acid concentrations were measured with a NanoDrop® spectrometer. RNA sequencing was performed using an Illumina Hiseq 2000 platform and TruSeq v3 kits. Four µg of nucleic acids were used for RNAseq and the rest of the material was used for the quality control. RNA libraries were sequenced on a single flow cell single stranded 51 bp read. Sequencing reads have been submitted to the European Nucleotide Archive (ENA) and are available under the submission number PRJEB15282. The amount and the length of reads obtained from each analysis are summarized in the Supplementary Table 2.

Relative expression estimation. Gene copy numbers and expression levels were estimated from DNA and RNA read counts respectively, mapping to annotated ORFs using reference gene annotation from TriTrypDBv7. The actual read-counts were measured with BedTools (v2.19) [START_REF] Quinlan | BEDTools: a flexible suite of utilities for comparing genomic features[END_REF] using the BAM files obtained after mapping. In order to normalize for gene length, read counts were estimated in Reads Per Kilobase Mapped (RPKM) [START_REF] Mortazavi | Mapping and quantifying mammalian transcriptomes by RNA-Seq[END_REF] across all genes and samples. These RPKM measurements were then used to compare expression levels across samples. show apparent regimens of mutation rates. The green line fits a linear model for chromosomes with an aneuploidy frequency higher than 25%, while the blue line indicates chromosomes that are lower in frequency. Chromosome

Figures

  using the following formula With t isolate trisomic for chromosome N -T all the isolates trisomic for chromosome N d isolate disomic for chromosome N -D all the isolates disomic for chromosome N -MA(s,T) the allele s that has the highest average frequency across all the T isolates trisomic for Chromosome N -AvgFrequency(MA(s,T),T) is the average frequency of MA(s,T) across all the T isolates trisomic for Chromosome N -AvgFrequency(MA(s,T),D) is the average frequency of MA(s,T) across all the D isolates disomic for Chromosome N This measure reflects the tendency of a given allele to dominate in frequency the overall trisomic population. It makes it possible to discriminate between a scenario where the two original disomic alleles are equally likely to have increased frequencies in each isolate (i.e. random selection of the duplicated homologous chromosome), and situations where the same allele is consistently amplified across all isolates (i.e. clonal amplification of the same trisomy or consistent selection of the same homologous chromosome amplification). The comparison between SV and LOH across field isolates is shown in Fig1d, with a p-value calculated as the t-test statistic based on Pearson's product moment correlation coefficient using cor.test() function in R.

Figure 1 :

 1 Figure 1: Genome instability, aneuploidy co-occurrence and haplotype selection in 204 L. donovani field isolates. Polysomy analysis based on read-depth comparison and derived allele frequency analysis on specific aneuploidies is shown. (a) Polysomy analysis. Polysomy level was estimated for each chromosome and sample by read-depth analysis. The distribution of chromosome copy number is shown across all field isolates. (b) Cooccurrence analysis. Pearson correlation measured for each pair of chromosomes was plotted while considering the median read depth measured in every field isolate. Empty positions correspond to either uncorrelated pairs or insignificant correlations (significance level of 0.05). One main cluster appears containing typically trisomic chromosomes. (c) Allele frequency profiles. Merged allele frequency distributions are shown for chromosomes 1 and 5 for both disomic (right) and trisomic (left) states. (d) Segregation variation vs allele diversity analysis. This panel summarizes the relationship between segregation variations (SV, reflecting bi-modality of the merged frequency profiles) and loss of heterozygosity (LOH, reflecting allelic diversity) when comparing trisomic and disomic isolates across 204 field isolates and the 36 chromosomes of the genome. For most chromosomes, both quantities are highly correlated as indicated by the blue regression line.

Figure 2 :

 2 Figure 2: In vivo aneuploidy dynamics. The three panels document the dynamic nature and the mosaicism of L. donovani aneuploidy in situ during hamster infection. (a) Aneuploidy pattern of L. donovani field isolate BPK282/0 after over 20 passages in culture (in vitro, iv)

Figure 3 .

 3 Figure 3. In vitro aneuploidy dynamics. This figure shows read-depth variations in L. donovani LD1S genes (8760 genes) across chromosomes at different passages during culture adaptation (4 samples) and in the 8 sub-clones derived from the parasite population at passage 20. Chromosome gene read-depth distributions are shown in boxplots depicting upper and lower quartiles and the median. (a) Read depth variation during culture adaptation. Variations in read-depth for L. donovani LD1S amastigotes isolated from infected hamster spleen (sp-ama) and derived promastigotes at passage 2, 10, and 20 (p2, p10, p20). Only aneuploid chromosomes are shown (see full panel in Supplementary Fig. 3a). Normalized read depth for every gene within each chromosome is displayed using a standard box-plot representation with the central line representing the median of the distribution. (b) Read depth variation on individual sub-clones. Variations on these same

Figure 4 :

 4 Figure 4: Fluctuations of allele frequency during culture adaptation. These panels document haplotype selection by following the allele frequency distribution in 8 sub-clones and the original p20 parasite population (n=9 samples in total) across . (a) Haplotype selection in clones. Variable sites in chromosomes undergoing amplification were either colored according to the major alleles (i.e. highest frequency) or left in gray for balanced heterozygote sites. Each painted pattern represents a chromosome, either from the original p20 population (first line) or one of the 8 derived sub-clones. Allele frequency distributions for each sample are shown on the left of each panel. The color code of the distributions corresponds to the read depth, with disomic-to-trisomic-to-tetrasomic transitions indicated by a continuum colored from red to green. The number of heterozygous sites for each chromosome and sample are shown by the histograms (right side of each panel). Samples with low counts as observed for chromosomes 15 and 20 represent loss of heterozygosity occurring during culture adaptation. (b) Haplotype selection during culture adaptation. Allele frequency distributions corresponding to hamster-derived splenic amastigotes (sp-ama) and derived promastigotes during culture adaptation at passage p2, p10 and p20 are shown for selected chromosomes. All chromosomes undergoing polysomy variations are displayed in Supplementary Fig. 3.

Figure 5 :

 5 Figure 5: Aneuploidies and haplotype selection correlate with phenotype and fitness. This figure shows the relationship between aneuploidy and haplotype selection with

Figure 6 :

 6 Figure 6: Aneuploidy influence on heterozygosity. Each chromosome is represented by its number on the graph and colored according to aneuploidy frequency after culture adaptation, as measured in the 204 L. donovani field isolates. The two lines (blue and green)

  31 is excluded from the fitted

		10	21	17	24	25	36	19	34	31	32	27	30	9	20	5	26	11	14	22	8	23	1	6	3	28	12	13	15	2	16	7	33	18	29	4	35
	10	o																																			1
	21		o
	17			o
	24				o																																0.8
	25					o
	36						o
	19							o
	34								o																												0.6
	31									o
	32										o
	27 30											o	o																								0.4
	9													o
	20														o
	5															o																					0.2
	26																o
	11																	o
	14 22																		o	o																	0
	8																				o
	23																					o
	1																						o														-0.2
	6																							o
	3																								o
	28																									o
	12																										o										-0.4
	13																											o
	15																												o
	2 16																													o	o						-0.6
	7																															o
	33																																o
	18																																	o			-0.8
	29																																		o
	4																																			o
	35																																				o	-1

Data availability. All Sequencing reads have been submitted to the European Nucleotide Archive (ENA) and are available under the Accession Number PRJEB15282. The amount and the length of reads obtained from each analysis are summarized in the Supplementary Table 2.

Code availability

Custom scripts and code necessary to generate figures and files can be found in the following repository: https://github.com/pprieto/genomevoleish

HTSeq analysis to correlate read-depth levels between genome and transcriptome (Fig. 5a, Supplementary Fig. 7) and found a very high correlation for most chromosomes (r=0.72), with the notable exception of chr 31. In this chromosome that is predominantly tetrasomic across various Leishmania strains, transcriptomic output appears halved, suggesting some form of dosage compensation, which cancels the polysomic effect. We also found karyotypic fluctuations resulting from culture adaptation to correlate with rapidly increasing in vitro fitness as judged by the decreasing generation time (Fig. 5b, left panel). Concurrently with these variations the parasites showed decreased in vivo fitness as indicated by lower infectivity (Fig. 5b middle and right panels). To assess if these variations take place during the parasite life cycle, we separately sequence amastigotes obtained from spleen and liver from one infected hamster. We found significantly distinct allelic profile variations for chr 20, which showed a unimodal profile in liver and a bimodal profile in spleen (Fig. 5c, Supplementary Fig. 8), while read depth analysis indicated a disomic state in both tissues (Supplementary Fig. 9). These observations closely recapitulate our in vitro analyses of parasite clones and indicate that transient polysomies may contribute to parasite adaptation in vivo.

Long-term genetic implications on genetic diversity

Maintaining diversity under strong selection poses a dilemma for all microbial pathogens, since it requires a mechanism that permits some compromise between immediate survival through selection and long-term adaptation through the maintenance of genetic diversity. We therefore searched for traces of increased genetic diversity associated with frequent aneuploidies and found a clear signal in the 204 field isolates (Fig. 6). Integrating genetic

Authors

Competing interests

The authors declare no competing financial interests. models as it stands out as the only chromosome with aneuploidy in all samples. 765