A complete Leishmania donovani reference genome identifies novel genetic variations associated with virulence

Leishmania donovani is responsible for visceral leishmaniasis, a neglected and lethal parasitic disease with limited treatment options and no vaccine. The study of L. donovani has been hindered by the lack of a high-quality reference genome and this can impact experimental outcomes including the identification of virulence genes, drug targets and vaccine development. We therefore generated a complete genome assembly by deep sequencing using a combination of second generation (Illumina) and third generation (PacBio) sequencing technologies. Compared to the current L. donovani assembly, the genome assembly reported within resulted in the closure over 2,000 gaps, the extension of several chromosomes up to telomeric repeats and the re-annotation of close to 15% of protein coding genes and the annotation of hundreds of non-coding RNA genes. It was possible to correctly assemble the highly repetitive A2 and Amastin virulence gene clusters. A comparative sequence analysis using the improved reference genome confirmed 70 published and identified 15 novel genomic differences between closely related visceral and atypical cutaneous disease-causing L. donovani strains providing a more complete map of genes associated with virulence and visceral organ tropism. Bioinformatic tools including protein variation effect analyzer and basic local alignment search tool were used to prioritize a list of potential virulence genes based on mutation severity, gene conservation and function. This complete genome assembly and novel information on virulence factors will support the identification of new drug targets and the development of a vaccine for L. donovani.

Second-generation sequencing technologies including Illumina, have made the sequencing of large genomes feasible through the mapping of short sequence reads of 50 to 250 nucleotides (nt) to a reference genome 5 . While human and many other well studied higher vertebrates have better assembled reference genomes 6 , the kinetoplastids suffer in this regard because most Leishmania species either lack sequencing information altogether or have incomplete reference genomes with sometimes thousands of sequence gaps 7 . The current L. donovani reference genome (ASM22713v2 from strain BPK282A 8 ) was generated using second generation technologies and contains over 2,000 gaps and therefore there are many incomplete or inaccurate protein coding sequences. The first complete Leishmania genome generated is that of L. major by a consortium of laboratories employing large insert clone tiling paths to sequence each chromosome individually 9,10 . This genome was later improved by the reassembly of complex collapsed loci that were incorrect in the original reference genome 11 .
Since then however, advances in sequencing technologies have drastically reduced the cost of sequencing and eased genome assembly tasks by increasing the length of the generated sequences. Long read sequencing or "third-generation" sequencing refers to more recent technologies including Oxford NanoPore 12 and Pacific Biosciences (PacBio) 13 that can result in reads ranging up to 50 kb or 100 kb that are capable of generating more complete genomic assemblies, provided the read lengths traverse across repetitive elements. One such highly repetitive cluster is the A2 gene family from L. donovani considered to be an important virulence factor and is necessary for survival in visceral organs 14,15 and protection against host response stress 16,17 . Due to its repetitive nature, the A2 gene cluster is misassembled in all Leishmania genomes generated using second generation sequencing, and only resolved in a recent resequencing effort targeted to L. infantum exploiting the long-read capabilities of PacBio sequencing which resulted in a complete genome assembly 18 . The current L. donovani genome however was obtained from second generation sequencing and consequently, no precise DNA or complete protein sequences are available for any A2 protein in L. donovani, hindering the comparison of A2 genes in visceral disease-causing strains or using mass spectrometry to identify A2 proteins which relies on accurate genome sequences for protein identification.
In this study, we have combined second and third generation sequencing to generate a complete assembly of the L. donovani genome from the strain responsible for cutaneous leishmaniasis (CL) in Sri Lanka 4,19 . This new assembly enabled the generation of an improved genome annotation and an unbiased analysis of chromosome synteny comparing L. donovani and L. major genes and strand switch transcription units. We have used this complete assembly to re-interrogate the genetic makeup of the visceral and cutaneous disease-causing L. donovani strains resulting in the identification of novel SNPs and indels generating a more complete and accurate chromosome map of the genetic differences between these phenotypically distinct L. donovani strains 4,20 . This study further enabled re-annotation of much of the genome highlighting the importance of a complete reference assembly to support future functional genomic and proteomic studies involving the L. donovani pathogen.

Results
A complete L. donovani genome assembly. The currently available assembly for L. donovani (ASM22713v2 from strain BPK282A 8 ) contains over 2,000 gaps due to the presence of low complexity regions and the highly repetitive nature of the Leishmania genome 21 . This incomplete assembly makes it difficult to compare L. donovani genomes from strains with different phenotypic properties. DNA was therefore isolated from the attenuated cutaneous disease-causing strain of L. donovani from Sri Lanka 4 and was subjected to deep sequencing using second and third generation sequencing. We reasoned that a complete assembly of the genome from this attenuated L. donovani strain will identify a more complete complement of genetic changes associated with loss of virulence of this strain. A total of 9 PacBio sequencing runs were performed generating 712,443 reads representing an estimated 107-fold coverage of the estimated 35 Mb genome. Importantly, there were 51,484 reads longer than 12 kb, representing a 20-fold coverage in very long reads. The long-read sequencing data was assembled using various assemblers as described in methods and merged using the longest chromosomes produced by each assembler followed by refinement using the high-quality short-read Illumina-generated data and iterative edge extension to close the remaining gaps.
The previous L. donovani reference assembly (ASM22713v2 from strain BPK282A) had over 2000 gaps spread across the 36 chromosomes. Figure 1 depicts the location on each chromosome of the gaps that have been closed in the new assembly reported here. The new assembly now contains contiguous DNA sequences in all 36 chromosomes and a corresponding 22-fold increase in N50 indicating that a larger proportion of the data has been assembled into large contigs as 50% of the genome is contained in contigs >= N50, resulting in an N50 of over 1Mbp (Table 1). Further, using this completed assembly, we have generated annotations for more potential protein coding regions than previously annotated (8,633 compared to 7,969 proteins) and identified more transfer-RNA and ribosomal RNA genes as well as all 6 small nuclear RNA genes, all spliced leader RNA genes and close to a thousand small nucleolar RNA genes. An additional 13 genes were marked as pseudogenes due to multiple internal stop codons and/or frameshifts. (Supplementary Table S1). Alignment of the second-generation Illumina reads to the PacBio generated assembly was used to cross-validate and correct the assembly at the nucleotide level. Graphs of the coverage from the alignment of Illumina and PacBio data across the 36 chromosomes are available in Supplementary Fig. S1. Taken together, we consider this new assembly to be contiguous and complete.
Assembly of the A2 virulence gene cluster and synteny comparison between L. major and L. donovani. A2 is a major virulence factor required for L. donovani survival in visceral organs 22 . The A2 gene family cluster on chromosome 22 has recently been assembled for L. infantum 18 , however has not been for L. donovani. We therefore investigated whether the structure of this region could be determined with this revised assembly. It was advantageous that the attenuated cutaneous L. donovani strain used in this assembly has fewer copies of the A2 genes than other virulent strains of L. donovani 4 . As shown in Fig. 2a, the new assembly could read-through the entire cluster of highly repetitive A2 and flanking sequences and could position the A2 genes and interspersed flanking 3′ A2-rel and 5′ A2-rel genes. The A2 genes are contained in two opposite facing clusters on either side of a strand-switch locus consisting of one cluster of 3 copies and one cluster with a single A2 gene.
The long sequence reads generated by the PacBio sequencing were crucial in generating the assembly of the A2 genes where reads of 11 kb and longer are shown spanning the repetitive cluster (Fig. 2b).
To generate supporting evidence for this A2 gene assembly, Western blot analysis of the A2 proteins from this strain was performed to compare the number and sizes of the A2 proteins with the predicted molecular weights from this assembly (ORFs; Supplementary Fig. S2). As shown in Fig. 2c, the apparent molecular weights from Western blotting correspond to the sizes predicted from the sequenced ORFs. The 3 bands on the Western blot are consistent with the molecular weights of the 4 gene products as the A2 gene copies 2 and 3 encode proteins of a similar size ( Supplementary Fig. S2). This represents the first complete structure and sequence for A2 genes in L. donovani, a prototype virulence factor. The difficulty in assembling this complex region is demonstrated in Fig. 2d, where a deviation from the average read coverage can be seen around the 300,000 bp position, in and around the A2 cluster, due to difficulties in the aligner assigning a unique position to similar reads across a repetitive region.
Directly comparing synteny at the chromosomal level was not possible with the previous L. donovani assembly due to the heavy fragmentation of the genome. With the new L. donovani assembly, it was possible to accurately compare chromosome synteny between L. donovani and L. major. As shown in Fig. 3, the genome of L. donovani, exhibited a very high level of synteny with the L. major. Chromosome 22 was highlighted here because this is the location of the A2 genes that have become pseudogenes in L. major and have therefore diverged between these old-world species 15 . The level of synteny demonstrated here for chromosome 22 was maintained on all other chromosomes ( Supplementary Fig. S3). These results indicate that evolution between cutaneous and visceral pathologies by different Leishmania species resulted largely from SNPs, pseudogenes and copy number variation and not from large changes such as chromosome rearrangements or complete gene deletions/insertions.

Identification of new genes and improvements in annotations.
As this assembly was larger in terms of total number of bases covered and more contiguous due to the removal of sequence gaps, the impact this had on gene annotations was investigated. The genome from the new assembly was annotated using the Companion pipeline 23 and the new and previous annotations (GenBank: GCF_000227135.1) were then aligned together and overlapping annotations were removed. Remarkably, close to 15% of the L. donovani protein coding genes had new or edited annotations as shown in Fig. 4a. Part of this increase in number of annotations resulted from the expansions of multi-copy gene families beyond the copy numbers in the previous annotation. An example is shown in Fig. 4b where there are 10 amastin genes identified in this new assembly compared to the previously identified 2. These results support the use of this assembly as the reference for bioinformatic analysis as it provides a more complete and accurate annotation of the L. donovani genome ORFs.  Comparison of virulent and attenuated L. donovani parasites. As indicated above, there are 2 distinct strains of L. donovani in Sri Lanka where one is responsible for visceral leishmaniasis (VL) and the other for cutaneous leishmaniasis (CL) 4 . Subsequently, the CL strain was experimentally passaged through the visceral  organs of BALB/c mice to select for a gain-of-function strain with increased virulence (IV strain) for survival in visceral organs where it was revealed through proteomic analysis that the resulting IV strain had an increase in stress response and antioxidant proteins 24 . Illumina whole genome sequencing and comparative genomic analysis of the VL, CL and IV strains was performed to identify SNPs associated with a change in virulence for survival in the visceral organs. As shown in Fig. 5, all 70 of the previously identified homozygous SNP differences between the VL and CL strains 4 were confirmed in this new assembly and an additional 15 novel SNPs within protein coding genes were found using this complete assembly. In addition, there were 12 mutations associated with the IV strain with gain-of-function for increased survival in visceral organs not labeled in Fig. 5; four were heterozygous but with frequencies changing towards the VL genotype (IV → VL), four were heterozygous but present only in the IV strain and four were homozygous deletions in the IV strain. The newly identified differences between the VL and CL strains and the ones contributed from the IV strain are summarized in Table 2.
Combining data from the previous and current analysis, all the genes with genetic differences were organized into priority clusters based on the likelihood to affect protein function and phenotype (Fig. 6). A detailed list of the genes and cluster assignments is shown in Table 3. From the 66 genes containing 70 SNPs previously identified, 7 were previously experimentally assessed using gene replacement with a wildtype copy for virulence in visceral organs 4 and one gene was identified as a misannotation and was therefore removed from the list. In decreasing order of priority, 13 genes in the highest impact cluster (red) were characterized as potentially having a major effect on protein function due to large amino acid changes or co-occurrence of mutations in both the VL and IV strains. SNPs in common between the IV and VL strains (4 IV → VL in the red cluster) indicate a selection associated with survival in visceral organs during the experimental passaging of the CL strain in mouse visceral organ 24 . Nine genes either with multiple co-occurring SNPs or non-conservative amino acid changes in conserved domains with a high score as assessed by PROVEAN software were placed in the second highest priority cluster. As detailed in methods, PROVEAN is a bioinformatic tool that classifies the significance of specific genetic mutations with respect to protein function 25 . Eighteen genes with non-conservative amino acid changes occurring in conserved domains but scored as unlikely to have a large effect on protein function by PROVEAN were placed in Cluster 3. Twenty-one genes with conservative amino acid changes in conserved domains were divided between Cluster 4 and Cluster 5 based on PROVEAN scores and 14 genes with mutations in domains with higher variability were placed in Cluster 6. Four mutations seen solely in the IV strain but not in the VL or CL strains are likely to be the result of random mutation or adaptation specific to the murine host were placed in Cluster 7.
A 25 kb region on chromosome 36 containing 4 genes was found to be missing in the IV strain but present in the VL and CL strains. This deletion did not occur in a location previously identified on this chromosome where a fission can occur as seen in L. alderi 26 . Upon experimental verification, this deleted region was present in a Compared to the previous L. donovani reference assembly (ASM22713v2 from strain BPK282A1), there were 1,087 protein coding genes unannotated or differently annotated in the current assembly. Unannotated or differently annotated genes were obtained by removing all annotations generated from our assembly that shared 95% or greater similarity to those previously available 8 . (b) Expansion of the amastin gene cluster on chromosome 8. Top track contains the previously two known coding sequences aligned to the previous L. donovani reference assembly (ASM22713v2 from strain BPK282A1). Gaps in the previous assembly depicted as dotted lines. Bottom track contains 10 amastin genes identified in the updated assembly. One previously identified Amastin gene has been aligned, 1 has been expanded and 8 have been annotated de novo.  Fig. S4). The enrichment of this deletion in the IV strain could therefore be a consequence of selection in the mouse and likely to be unrelated to human visceral disease because this region is present in wild type or VL strains of L. donovani as well as L. major and therefore classified in cluster 7.
The classification of genetic differences in the CL, VL and IV genomes summarized in Fig. 6 and Table 3 represents a prioritization of genes to be empirically investigated for controlling the different phenotypes of these virulent and attenuated strains.

Discussion
It has been possible to generate a complete genome assembly for L. donovani through combining second and third generation sequencing technologies, similarly to a recent resequencing of the L. infantum genome resulting in a complete assembly, highlighting the usefulness of PacBio sequencing in regards to Leishmania genomes 18 . This resulted in a more accurate annotation of the genome increasing the number of potential protein coding genes and identifying novel mutations/polymorphisms associated with virulence. It was remarkable that the present assembly resulted in annotation changes in close to 15% of the genome representing 1087 protein coding genes. Although 13 degenerate pseudogenes are identified in Supplementary Table S1 more do exist since our annotations derived from functional genes in L. major and therefore genes functional in other species were not identified. Through this updated genome annotation, additional SNPs have been identified including in genes potentially involved in visceral disease and several non-coding genes have been annotated allowing future L. donovani research beyond protein coding genes. It has also been possible to assemble known virulence factor gene families in L. donovani including the A2 and Amastin gene families. This version of the L. donovani genome assembly will significantly improve genomic, functional genomic and proteomic research outcomes and support the identification of drug targets and the development of vaccines. This assembly further provides a larger repertoire of target DNA sequences to identify diagnostic and prognostic disease progression markers. Given the recent interest in generating genetically modified live attenuated parasites as vaccine candidates 27 , a complete genome assembly will permit the verification that genetic modifications target intended genes with no off target mutations.
Supported by a de novo assembly, this study provides the first direct evidence for synteny between chromosomes in L. donovani and L. major, two old world parasite species with different pathologies and reservoirs. Previously, due to the large number of gaps in the L. donovani genome, the segments were aligned to a reference assembly assumed to be syntenic and only gene synteny was possible. In contrast, the contiguous assembly presented within used an entirely reference-free and by extension, bias-free generation process. This assembly can be used in future sequencing efforts aimed at comparing genes and synteny of genomes of other Leishmania species with L. donovani. The strong gene level synteny further highlights the major phenotypic effects of SNPs and indel mutations when comparing genomes from L. donovani strains causing visceral and cutaneous pathologies. As no major chromosomal rearrangements or deletions are apparent between phenotypically different Leishmania species as previously reported 2,4 , and including this study, suggesting that virulence and tropism can be acquired or lost through relatively small coding changes at the amino acid level such as SNPs, indels and frameshifts without the need for chromosomal scale events. This study reports the complete A2 gene continuous sequence and an assembly of an entire A2 cluster including the A2-rel flanking genes in L. donovani. While the organization of this A2 gene cluster was previously theorized based on available sequences and Southern blot analysis, no sequencing technology could accurately read-through an entire cluster 4,28 prior to third-generation sequencers. Interestingly, a similar organization of A2 and A2-rel flanking genes was obtained during the resequencing and assembly of the L. infantum genome 18 , further supporting the genomic arrangement of this important virulence cluster. The present assembly contains entire A2 ORFs that are consistent with the corresponding protein sizes determined by Western blot analysis and provides novel insight into this elusive virulence factor through identifying 2 amino acid insertions between the 10 amino acid repeats at geometric intervals as well as a defined C-terminal sequence (Supplementary Fig. S2). The deviations from the 10-amino acid repeat sequence could contribute to the proper folding and function of the A2 protein.
In an attempt to identify additional genes associated with survival in visceral organs, the attenuated cutaneous L. donovani strain was experimentally passaged continuously through the visceral organs of BALB/c mice over an 8 month period to generate a gain-of-function strain with increased survival in the liver and spleen and was termed the IV strain 24 . Sequence analysis of the IV strain in this study did not identify homozygous SNP differences with the parental cutaneous strain but did identify four heterozygous SNPs with the same sequence in the virulent L. donovani strains, classified as high impact in Fig. 6. The corresponding SNP-containing genes with unknown function are of high priority for future studies. Nevertheless, although the gain of function IV  Table 2. Summary of novel mutations identified in this study. All mutations are annotated using VL as the wild type amino acids and CL as the mutated amino acids. Genes with annotations in the previous assembly list the previous gene ID in italic, genes annotated only in this assembly list only one gene ID. The top segment lists fifteen attenuated cutaneous strain specific mutations identified in this study. Mutations marked with a appear at 50% but also co-occur with gene duplication event and are therefore possibly homozygous on one copy. strain had significantly increased survival in visceral organs 24 , it was surprising that this strain did not have more genetic differences associated with the increased virulence. It is possible that the selection process for survival in the visceral organs of mice is different from that in humans. The Illumina sequence analysis of the cutaneous (CL) and visceral (VL) disease associated strains using the complete assembly identified 15 novel homozygous SNPs beyond the previously identified 70 SNPs 4 (Fig. 5). One of these new SNPs was in the Raptor gene that is part of the highly conserved Target Of Rapamycin (TOR) signaling pathway 29 . There are three TOR gene homologs in the Leishmania genome 30 revealing this pathway is conserved in kinetoplastids. Interestingly, the RagC GTPase which is a binding partner of Raptor in the TOR pathway is also mutated in the attenuated cutaneous L. donovani strain and restoration of the wildtype RagC GTPase increased virulence in visceral organs 4 . Considering that there are two mutated genes (RagC and Raptor) in the TOR pathway in the attenuated cutaneous L. donovani strain strongly highlights this pathway as playing a role in determining disease tropism and virulence.
As both HIVE and VarScan were used to identify SNPs and indels, we are confident that the expanded list of 83 variable genes shown in Fig. 6 contains most if not all the genes associated with visceral disease, with the exception of UTR mutations that may influence protein expression levels. Since this number of genes is relatively small, we are currently investigating all genes in clusters 1-4 with respect to their involvement in visceral organ virulence using CRISPR-Cas9 gene editing recently developed for use in Leishmania 31,32 . It is noteworthy that the correct selection of gRNA sequences for CRISPR-Cas9 gene editing requires a complete genome and accurate annotations for precise gene editing with no off-target mutations that is now possible with the complete assembly reported here.

Methods
Whole genome sequencing. DNA extraction. Leishmania DNA for both Illumina and PacBio sequencing was derived from the attenuated cutaneous strain of L. donovani from Sri Lanka 4 that was passaged through mice to increase survival in visceral organs 24 . DNA was extracted following the previously described phenol-chloroform methods for isolation of Trypanosomatid genomic material 33 .
PacBio sequencing. A total of 9 sequencing cells were prepared. 7 cells were prepared using the DNA Template kit v2.0 (3-10 kb) with DNA/Polymerase Binding Kit P4 and 2 using the DNA Template Prep Kit 3.0 with DNA/ Polymerase Binding Kit P5. The cells were sequenced on a PacBio II RS instrument with BaseCaller v1 protocol. Genome assembly. HGAP assembler. Raw reads from the 9 sequencing cells were loaded into the SMRT Analysis portal (Pacific Biosciences) in HD5 format. The Hierarchical Genome Assembly Process (HGAP) version 2 with Quiver polishing was chosen as version 3 is stated to improve speed at the detriment of assembly quality. Expected genome size was set to 36Mbp, minimum read length for pre-assembly was set to 500 bp and minimum read length for full assembly was set to 100 bp. Minimum Polymerase Read quality was set to 0.80, and the remainder of options remained at default settings. Figure 6. Summary of all genes with non-synonymous mutations between the cutaneous, visceral, and gainof-function strains of L. donovani. All non-synonymous SNPs and Indels were classified as common to our previous study (2014 CL 4 ) or identified in this study (Novel), as well as by their effect on amino acid changes from top to bottom, colored red to green in descending order of likelihood to affect the phenotype of the parasite. 66 genes were common to the previous data set. Of those genes, 7 were previously investigated 4 and 1 was rejected due to an open reading frame misannotation. 25 genes were only listed in this study (Novel). Diagram created using SankeyMATIC (http://sankeymatic.com). LdCL_310020800 LdBPK_311390 LdCL_250013200 LdBPK_250790 LdCL_310020100 LdBPK_311320 LdCL_310022200 LdBPK_311510 LdCL_080011700 LdBPK_080670 LdCL_340029800 LdBPK_342210 LdCL_230014900 LdBPK_230830 LdCL_310037100 LdBPK_312870 LdCL_220015800 - Cluster 2 (9) Multiple SNPs in the same gene, Non-conservative amino acid change in conserved region with good PROVEAN score LdCL_270015000 LdBPK_270840 LdCL_290026900 LdBPK_292100 LdCL_310021600 LdBPK_311470 LdCL_310028800 LdBPK_312080 LdCL_340046300 LdBPK_343690 LdCL_290022800 LdBPK_291720 LdCL_310024300 LdBPK_311710 LdCL_360006000 LdBPK_360120 LdCL_360062000 LdBPK_365480 Cluster 3 (18) Non-conservative amino acid change in conserved region with poor PROVEAN score LdCL_070018300 LdBPK_071330 LdCL_320013800 LdBPK_320820 LdCL_250016900 LdBPK_251150 LdCL_220022000 LdBPK_221470 LdCL_250015300 LdBPK_251000 LdCL_040011100 LdBPK_040560 LdCL_360016300 LdBPK_361120 LdCL_200014300 LdBPK_200960 LdCL_090011700 LdBPK_090660 LdCL_130016200 LdBPK_131090 LdCL_340009000 LdBPK_340390 LdCL_130017800 LdBPK_131230 LdCL_230009900 LdBPK_230440 LdCL_220018100 LdBPK_221070 LdCL_340044900 LdBPK_343550 LdCL_290028400 -LdCL_250011400 LdBPK_250620 LdCL_270014900 LdBPK_270830 Cluster 4 (8) Conservative amino acid change in conserved region with good PROVEAN score LdCL_350013100 LdBPK_350830 LdCL_360052700 LdBPK_364550 LdCL_230026600 LdBPK_231940 LdCL_230009400 LdBPK_230400 LdCL_320031100 LdBPK_322560 LdCL_020008200 LdBPK_020280 LdCL_310027700 LdBPK_311990 LdCL_320031200 LdBPK_322570 Cluster 5 (13) Conservative amino acid change in conserved region with poor PROVEAN score LdCL_330011900 LdBPK_330640 LdCL_170010200 LdBPK_170470 LdCL_070011900 LdBPK_070700 LdCL_210025000 LdBPK_211930 LdCL_290022900 LdBPK_291730 LdCL_200006300 LdBPK_200140 LdCL_290029000 LdBPK_292290 LdCL_030007500 LdBPK_030250 LdCL_250006200 LdBPK_250110 LdCL_360015800 LdBPK_361070 LdCL_310023400 LdBPK_311630 LdCL_360062700 LdBPK_365540 Celera assembler. The PacBio corrected Reads (PBcR) module of the celera-assembler version 8.3 was used to assemble the long reads data 34 . The subreads were first extracted from the PacBio H5 files to FASTQ using bash5tools.py. The Bogart unitigger was used by specifying the "unitigger = bogart" option in the spec file. The consensus caller module was PBDAGCON. Due to the sequences originating from a non-clonal sample and the use of the DNA/Polymerase Binding kit P4 in some PacBio sequencing cells which produces lower quality data than P5 kits, error rate limits were relaxed for various variables, listed in the full spec file available in supplementary information (Supplementary Methods S1).
Canu assembler. The Canu v1.0 assembler is a modified version of the Celera Assembler designed to handle high noise data such as NanoPore and PacBio sequencing data. Canu has both the ability to assemble raw PacBio data by performing error correction using consensus sequence or assemble data in a hybrid mode where PacBio reads are pre-error corrected using short read Illumina data. In Raw mode, the trimmed PacBio reads were given to the assembler using default settings except for the expected genome size option which was set to 35Mbp using the option "genomeSize = 35 m". In hybrid mode, the Illumina reads were first error corrected by internal consensus using Pollux 35 , the paired end reads were then merged together to form longer sequences with a high confidence core region using FLASh 36 , and used to correct the PacBio reads using Proovread 37 . The error-corrected PacBio data was then used by Canu to generate a draft assembly.
Pilon. The Pilon error correcting software was used to fix small errors present in the PacBio based assemblies using high depth and high accuracy Illumina data 38 . The entire MiSeq dataset in FASTQ format from the corresponding sample was aligned to the draft assembly using the Burrows-Wheeler Aligner (bwa) to generate SAM alignment files. Samtools was then used to convert and sort the files to a binary usable format as described in the Samtools section. This alignment was then passed to the Java Pilon executable for correction of small indels, SNPs, gap filling and assembly of unmapped reads using the command "java -Jar pilon.jar-genome [new-assembly. fasta] -frags [alignment.bam] -fix all,novel".
GMCloser. GMcloser was used to merge the assemblies generated using the different assemblers 39 . Short read Illumina data was aligned to the contigs resulting from the different contigs from different assemblers with identical reads mapped to them were assumed to be part of the same chromosome. When a contig from one assembly encompassed a gap present in another assembly, the gap was filled with the missing information to generate a merged assembly with the least number of gaps. All the alignment and merging steps are handled internally to GMcloser using the command "gmcloser -t [assembly1.fasta] -q [assembly2.fasta] -r miseq_R1.fastq miseq_ r2.fastq -et".

Cluster Number
Cluster Mutation Type New annotation Equivalents (when available) Cluster 6 (14) Non-conservative amino acid change in less conserved region, Conservative amino acid change in less conserved region LdCL_340022100 LdBPK_341580 LdCL_060011600 LdBPK_060650 LdCL_210015400 LdBPK_211040 LdCL_050010900 LdBPK_050580 LdCL_230011600 LdBPK_230610 LdCL_250024100 LdBPK_251840 LdCL_070015100 LdBPK_071060 LdCL_250005300 LdBPK_250040 LdCL_310022100 LdBPK_311500 LdCL_200006800 LdBPK_200200 LdCL_250014400 LdBPK_250910 LdCL_230010400 LdBPK_230500 LdCL_120008300 LdBPK_120275 LdCL_290028100 LdBPK_292210 Cluster 7 (4) IV-only mutations LdCL_140010000 LdBPK_140470 LdCL_310035800 LdBPK_312770 LdCL_310036400 LdBPK_312810 LdCL_320046000 LdBPK_324000 LdCL_360021300 LdBPK_361580 LdCL_360021400 LdBPK_361590 LdCL_360021500 LdBPK_361600 LdCL_360021600 LdBPK_361610 Table 3. Summary of all genes containing mutations in the cutaneous isolates and classification into clusters. Entries were not repeated in multiple lists. Identified mutations were further classified into priority clusters for effect on protein function and future analysis for genes associated with survival in visceral organs. Mutations were prioritized by likelihood of contributing to visceral tissue tropism by severity of the coding change, accumulation of secondary mutations and conservation. Gene loci listed from the current assembly as well as previous ID numbers when available. IGV. The Broad Institute Integrative Genome Viewer 40,41 was used to perform quality control on assemblies and manually inspect fragments in order to close gaps. The Pilon tools was used with the "-fix novel" option which assembled short contigs from unmapped data. The fragments were then placed on the appropriate likely chromosomes based on gene annotations and submitted to another round of gap filling using Pilon and GMCloser to find reads supporting this placement or were removed if no reads supported the join.
Annotations. Companion. The Companion webtool (https://companion.sanger.ac.uk/) was used to annotate genes on the assembly contigs and refine the assembly 23 . The closest available reference organism was chosen (L. major) with the following options: contiguate pseudochromosomes, align reference proteins to target sequence, perform pseudogene detection, use RATT Species transfer type, and the L. donovani taxon ID. Additional L. donovani and L. major genes not automatically transferred were manually verified and appended if necessary. An additional 3 genes were manually added from a search of all ribosomal protein transcripts in trypanosomes. The snRNAs U1 through U6, ribosomal RNAs and the spliced leader RNA were manually annotated as necessary from the sequences available for L. major on TriTrypDB 42 . Sequences for H/ACA and C/D box snoRNA were manually mapped using published L. major snoRNA research 43 .
Galaxy. The Galaxy webtool (https://usegalaxy.org/) 44 was used to perform file conversions and data extraction such as moving a chromosome's FASTA sequence from one assembly to another.
Identification of new genes. Genomic annotations from the Companion Pipeline were downloaded in General Feature Format (GFF) and gene annotations were extracted using the Galaxy tool "Extract features" set to look for the "CDS" keyword in column #3 of the GFF file. Known coding regions from the reference L. donovani strain BPK282A1, assembly ASM22713v2 were downloaded from GenBank and aligned to our improved assembly in BED format. Bedtools intersect intervals through Galaxy 45,46 was used to identify annotations that were unique to our annotations or were not at least 95% covered previously using settings "-wa -f 0.95 -v -r".
Synteny. The online SynMap2 software 47 was used to generate the synteny dotplot across the entire genome using annotations from L. major and the annotations generated by Companion in this study. The chromosome to chromosome circular charts were generated by Companion as part of the annotation process.

Comparison of visceral (VL), cutaneous (CL) and increased virulence (IV) L. donovani
strains. BWA. The Burrows-Wheeler Aligner (BWA) was used to process the FASTQ Illumina sequencing files obtained from Genome Quebec. The maximal exact match algorithm was used in paired-end mode using the command "bwa mem" and providing the matched pair read files and reference sequence as arguments in order to generate a SAM format alignment file of the reads on the reference 48 .
Samtools. The samtools package was used for file manipulations and conversions 49 . The commands "samtools view -b" was used to convert the BWA generated SAM file to the binary alignment BAM format. The file was then sorted by alignment location for compatibility with downstream analysis software using "samtools sort -@ 30 -o [output.file]". The alignment files where then prepared for analysis using the mpileup modules which tabulates the base distribution at every position using the command "samtools mpileup -B -f [reference assembly] [strain specific position sorted BAM file] > [output.file]".
SnpEff. To filter the VCF files generated by VarScan to a list of non-synonymous SNP, we used the SnpEff software 51 . The oriented and annotated assembly was downloaded from the Companion tool as described above along with the gene annotation file in GFF format containing the names, locations and amino acid sequences of identified genes. This GFF file was used to build a SnpEff database using the SnpEff.jar command "build" with argument "-gff " after installing the genome and GFF file in the appropriate locations according to the software instructions. The SnpEff software was then used to annotate the 10 th column of the VCF file with mutation effect codes. All the mutations were then examined manually for accuracy using the Integrative Genomics Viewer (IGV) with all raw Illumina data loaded.
Classification. Non-synonymous mutations were clustered according to the mutation effect in order to prioritize further gene function studies. Each cluster was further broken down based on the mutation's PROVEAN score 25 . The PROVEAN software was designed to predict the magnitude of a mutation's impact on protein function. To generate PROVEAN scores, we retrieved homologous sequences from other Leishmania species and kinetoplasts and generated a multiple sequence alignment (MSA). The MSA was then passed to the PROVEAN software which scored each SNP based on the alignment. We used PROVEAN scores below a threshold of −2.5 as an indication a SNP is likely to affect protein function. Cluster assignments were as follows: