Resequencing of the Leishmania infantum (strain JPCM5) genome and de novo assembly into 36 contigs

Leishmania parasites are the causative of leishmaniasis, a group of potentially fatal human diseases. Control strategies for leishmaniasis can be enhanced by genome based investigations. The publication in 2005 of the Leishmania major genome sequence, and two years later the genomes for the species Leishmania braziliensis and Leishmania infantum were major milestones. Since then, the L. infantum genome, although highly fragmented and incomplete, has been used widely as the reference genome to address whole transcriptomics and proteomics studies. Here, we report the sequencing of the L. infantum genome by two NGS methodologies and, as a result, the complete genome assembly on 36 contigs (chromosomes). Regarding the present L. infantum genome-draft, 495 new genes have been annotated, a hundred have been corrected and 75 previous annotated genes have been discontinued. These changes are not only the result of an increase in the genome size, but a significant contribution derives from the existence of a large number of incorrectly assembled regions in current chromosomal scaffolds. Furthermore, an improved assembly of tandemly repeated genes has been obtained. All these analyses support that the de novo assembled L. infantum genome represents a robust assembly and should replace the currently available in the databases.

SCIENTIfIC RepoRtS | (2017) 7:18050 | DOI: 10.1038/s41598-017-18374-y sequencing technologies 9 , together with a significant reduction of sequencing costs, has speeded up the development of sequencing projects. As a positive consequence of this, an exponential increase in the number and diversity of sequenced genomes is taking place. In particular, based on the use of these new technologies, genomic drafts for several Leishmania species (and strains) have been published [10][11][12][13][14][15][16][17][18] and/or are publically available in databases (e.g., TriTrypDB.org). Nevertheless, although these genomic sequences are providing valuable data, the L. major (Friedlin strain) genome, decoded ten years ago using classical Sanger sequencing, continues to be the best assembled genome for the genus Leishmania. In fact, given the remarkable degree of synteny observed between the genomes of the different Leishmania species, the L. major genome is being used as the reference for building the chromosomal scaffolds of the other Leishmania species. However, apart from the L. major (Friedlin) genome, the rest of Leishmania genomes sequenced to date must be considered as draft assemblies, taking into account that they have been assembled in a number of contigs (between 562 and 10305) that vastly outnumbers the haploid set of chromosomes (34)(35)(36), depending on the Leishmania species 19 ). The existence of a large number of repetitive DNA sequences, which are scattered along the Leishmania genomes [20][21][22] , prevents the complete assembly of the Leishmania chromosomes when using short-read sequencing approaches. Another challenging issue to resolve during the genomic assemblies is the precise determination of the gene copy number in loci consisting of multiple tandemly arranged identical genes, a common feature of the Leishmania gene organization 23 . In fact, as we recently reported, these issues were the cause of assembly collapses affecting seven genomic regions that were missed at the time of the L. major (Friedlin) genome assembly 24 .
For sequencing the genome of the L. infantum species, the JPCM5 strain (MCAN/ES/98/LLM-724) was selected 8 . This strain was first isolated in Madrid (Spain) from a naturally infected dog that developed VL 25 . This strain was cryopreserved at the Centro Nacional de Microbiología (Instituto de Salud Carlos III, Madrid, Spain) and distributed to different laboratories around the world. This strain showed a high degree of virulence when assayed for experimental infections of dogs 26 . Furthermore, it was found to be suitable for genetic manipulations, and importantly that the laboratory-modified parasites retained their virulence 27 . All these studies supported the selection of this strain for the generation of the first genome draft of a viscerotropic Leishmania species 8 . Genome sequences were produced by the whole-genome shotgun method, achieving a mean coverage of five-fold. The sequences reads were assembled into 562 contigs, which finally were grouped into 36 chromosomal scaffolds using as reference the L. major (Friedlin) genome 7 . At that time, the authors estimated that the missing genomic regions would account for at least 150,519 bases 8 .
Here, due to the clinical relevance of L. infantum together with the fact that this species is a widely used model for molecular biology purposes, we undertook the aim to improve its genome assembly, currently formed by several hundreds of contigs. For that objective, we sequenced the L. infantum (JPCM5 strain) genome using two different platforms, the Pacific Biosciences (PacBio) technology, used to produce sequencing reads of 10-15 kb in lengths, and the Illumina technology to generate paired-end short-reads. As a result, we are providing one of the most comprehensive reference genomes available, having a quality at least comparable to the L. major (Friedlin) genome assembly, and clearly better than the current L. infantum (JPCM5) one.

Results and Discussion
Illumina sequencing and assembly results. Total DNA isolated from L. infantum (JPCM5) promastigotes was sequenced with the Illumina high-throughput sequencing technology. A total of 56,327,604 paired-end 126 bp sequence reads with an average insert size of 310 bp were obtained. Thus, according to the estimated genome size of this Leishmania strain (32,134,935 bp 8 ), these reads would correspond to a sequencing coverage above 400×. A de novo assembly from these data was tried using several assemblers, and the best results were achieved with the CLC Genomics Workbench software (CLC Bio; version 5.0). After several refinements, the final assembly yielded 1874 scaffolds with a size of the longest contig of 363,515-bp and a total genome size of 31,179,733-bp. Thus, in spite of having a better sequence-depth, the assembly did not substantially improve the current draft for the L. infantum genome 8 . L. infantum and other Leishmania species contain many repetitive sequences, with sizes around 500-600 bp, scattered throughout the genome [20][21][22] . The high sequence identity shared by some of these repeated elements is a cause of conflicts for assemblers 28 when using short reads as those generated by Illumina and other sequencing platforms 9 . A strategy often followed to circumvent this drawback consists in ordering the contigs into scaffolds using a well-assembled genome from a related organism as the reference. In this context, the L. major (Friendly) genome has been used to assemble most of the Leishmania genomes currently available at the TriTryDB database. However, this is a risky strategy that could induce misassembling in divergent regions. In addition, a large number of tandemly repeated, multi-copy genes exist in the Leishmania genomes 23 , and to determine the exact copy number of those genes is another challenge to resolve, when working with short reads, since repetitive sequence regions tend to collapse into a single copy if no nucleotide differences exists to distinguish among the gene copies 28 . On the other hand, to define the copy number by comparison with a reference, assuming that the copy number is conserved among different Leishmania species, is not accurate.
PacBio sequencing is a key technology for Leishmania genome assembly. As short-read paired-end Illumina technology resulted insufficient to achieve a full assembly of the L. infantum genome, it was decided to generate long-read sequences based on the single-molecule real-time (SMRT) sequencing technology developed by Pacific Biosciences (PacBio 29 ). A total of 311,471 reads with an average length of 11700-bp were obtained. A coverage of around 100× was estimated, taking into account the size of the L. infantum genome currently available at GeneDB. Interestingly, de novo genome assembly using the RS_HGAP_Assembly.3 protocol of the Pacific Biosciences SMRT Analysis Software v2.3.0 30 yielded 85 contigs, a number closer to the real number of L. infantum chromosomes, i.e. 36 19 . Of these, 41 were discarded as "spurious or artefactual", some of which corresponded to maxicircle sequences and others were short sequences showing low read coverage and high sequence similarity with regions present in the long contigs (see Materials and Methods for further details).
SCIENTIfIC RepoRtS | (2017) 7:18050 | DOI:10.1038/s41598-017-18374-y From the 44 remaining contigs, with the aid of the contigs generated from the Illumina reads data (see above) and using different assemblers and other bioinformatics tools, it was possible to join some of them, and finally reducing the total number of contigs to 36. Chromosomes 7, 12, 15, 22, 26, 28, 33 and 35 resulted by joining two contigs, but the rest of chromosomes were directly assembled from PacBio reads as a sole contig. Figure 1 shows the pipeline followed to achieve the final assembly (further technical details have been provided in the Materials and Methods section).
The accuracy of the assembled 36 contigs was assessed by alignment of both Illumina and PacBio reads to the assembled chromosomes. The read distribution on those chromosomes formed by two PacBio-derived contigs is shown in Fig. 2. It was observed a homogeneous distribution of reads along the joined regions, suggesting that the assembly was properly performed.
In summary, our results demonstrated that PacBio sequencing is appropriate for achieving an effective assembly of Leishmania genomes, but Illumina sequencing was relevant for accurately joining some contigs and for extending the chromosomal ends of chromosomes ( Fig. 1). Improvements in the L. infantum genome assembly. Table 1 summarizes major changes introduced in the new genome assembly regarding the previous one 8 . The total size of the genome (32,802,969-bp) has increased in 680,199-bp, and the number of undetermined nucleotides has been reduced from 20,399 in the current genome to zero in the new assembly. In addition, the chromosome 0, a chromosome created by the artificial joining of 34 genomic regions that could not be assigned with certainty to any of the 36 chromosomes, has disappeared after the de novo assembly reported here.
In agreement with the increase in the genomic size, the number of annotated genes has increased in 601, growing from 8195 (ref. 8 ) to 8796 (including also non-coding RNAs), in the de novo assembly reported here. In particular, we have identified 495 new protein-coding genes ( Table 2). An important fraction of these new genes corresponds to tandemly repeated genes that collapsed into one or two copies in the assembly of the current reference genome due to the relatively short size of the cloned fragments (4-kb or lower) and the sequence lengths (600-800 bp 8 ). Thus, a significant increase in the number of annotated genes has been determined in the loci coding for cysteine peptidase B (CPB; LinJ.08), ATG8 (LinJ.09 and LinJ.19), GP63 -leishmanolysin (LinJ.10), alpha tubulin (LinJ.13), elongation factor 1-alpha (LinJ.17), glycerol uptake protein (LinJ. 19   3 and 4 illustrate the corrections introduced regarding the copy number in the loci for beta-tubulin and glucose transporter 2, respectively. For comparison, it is shown the structure of these loci in L. major (Friedlin) and in the current reference L. infantum (GeneDB.org) genomes, together with the one derived from the new assembly described here (panels A). The homogeneous distribution of reads (either from Illumina or PacBio data) along the new assembled genome (panels C) contrast with the irregular distribution found when the reference genome (GeneDB.org) was used for reads-alignment (panels D). In the current L. infantum assembly (LinJ-Ref), a sole beta-tubulin gene was assembled, whereas six genes were assembled in the L. infantum genome reported in this work (LinJ-New). Noticeably, 16 beta-tubulin genes are assembled in the L. major (Friendly) genome (see Fig. 3). By contrast, for the glucose transporter 2 the current (LinJ-Ref) and L. major genomes show three copies, but in the new assembly (LinJ-New) six copies were detected (see Fig. 4). For most of the loci containing tandemly repeated genes, it is likely that the correct number of genes has been assembled, as reads coverage was found to be quite homogeneous along the 36 chromosomes of the new assembled L. infantum genome (Fig. 2, and data not shown). The sole exception is the rDNA locus, for which our assembly contains three tandemly linked copies of the 18 S rRNA-5.8 S rRNA-24S α/β rRNA unit; however, by measuring the mean read depth of this genomic region and normalization with the mean read coverage of chromosome 27, it was determined that the total number of rDNA units would be between 6 and 9 copies. The size of repetition unit (~10.5 Kbp) justifies that the PacBio reads also collapsed during the assembly of the rDNA locus. In current L. infantum   Nomenclature of the new assembled genome and synteny analysis. Among the new genes identified in the de novo assembled genome (Table 2) there are some genes that were not annotated in current databases because the corresponding genomic regions were missed. In addition, the new genome sequence has allowed to complete partial gene sequences and correct some sequence uncertainties (Ns) or errors. Current ID nomenclature of L. infantum genes (GeneDB.org) has been maintained as much as possible, as it is widely used and there are many publications having relevant data referring to that nomenclature. Hence, the former ID names have been kept even when some genes were found to be located at a different chromosome in the new assembly, regarding current reference genome. When new genes were needed to be named, the same nomenclature rules were used, and intercalated ID numbers were assigned to name these new genes. Finally, a total of 75 previously annotated genes have been excluded from the annotation of the new L. infantum genome. Some of the eliminated ID names corresponded to genes that were annotated as tandem gene duplications, but such gene duplications were not found in the new assembly. Other previous annotations were based on genomic regions that were bound  in an artefactual manner (see below). In summary, we have tried to maintain as much as possible the former gene IDs, even when the ORFs had to be partially corrected. For a few cases, the former IDs were eliminated to avoid confusion, as the new ORFs were very different to the previously annotated. On the other hand, a comparison between the current genome and the de novo annotated one highlighted important reorganizations in most of the chromosomes. In Fig. 5, two illustrative examples, affecting chromosomes 7 and 13, are shown. In these cases, the synteny of the homologous chromosomes in the L. major (Friedlin) genome is also shown. Remarkably, the de novo assembled genome is more syntenic when compared to the L. major assembly than when compared to the L. infantum reference genome, further supporting the view that some of these segments were incorrectly placed during the current L. infantum genome assembly 8 . In summary, all these analyses support the conclusion that the de novo assembled L. infantum genome described in this work represents a robust assembly. Full synteny maps comparing L. major (Friedlin) genome and the L. infantum genomes (new and reference) are provided in Supplementary Figures S1-S36. Conclusions. The availability of a robust genome sequence is a valuable resource for studies addressing whole-organism aspects following either genomics, transcriptomics or proteomics approaches. Advances in sequencing technologies have greatly facilitated genome sequencing tasks. However, genome sequencing data by themselves may have a limited utility unless they are adequately assembled in order to define the arrangement of genes and genome architecture. An enormous effort was invested to elucidate the genome structure and sequence of L. major, a milestone achieved in 2005. Genome comparison among the Leishmania species has shown a high degree of genome conservation in terms of both gene content and gene synteny across the genus. This finding led to the use of L. major (Friedlin) genome as a template to facilitate the assembly process of most of the Leishmania genomes reported afterwards. However, as demonstrated in this and other works 32 , that approach may lead to introduce assembly errors that would compromise future studies regarding gene content, gene models and genome architecture. Here, we present a de novo assembly of the L. infantum (JPCM5) genome based on sequence data derived from both long (PacBio) and short (Illumina) reads that yielded the expected 36 chromosomal-size contigs, without discontinuities and undetermined sequence (Ns), which are abundant in the current genome (GeneDB.org). Furthermore, this work is providing a methodological pipeline to obtain a full closed genome of a Leishmania species (or related kinetoplastids).
The new L. infantum genome sequence and annotations will be available at EBI databases, and also at the Leish-ESP web site (https://leishseq.neocities.org/). Moreover, the complete annotation of the new genome is

Methods
Parasites and DNA isolation. The L. infantum reference strain, JPCM5 (MCAN/ES/98/LLM-724), was isolated by Dr J. Moreno's group (WHO Collaborating Centre for Leishmaniasis, Centro Nacional de Microbiología, Instituto de Salud Carlos IIII, Madrid, Spain) from a dog suffering from visceral leishmaniasis. Promastigotes were cultured at 26 °C in RPMI 1640 medium (Sigma-Aldrich), supplemented with 20% heat-inactivated foetal calf serum (Sigma-Aldrich). Genomic DNA was isolated following the classical phenol-chloroform-isoamyl alcohol extraction method as described previously 33 . Illumina sequencing and reads assembly. Library construction and paired-end library sequencing were performed at the Centro Nacional de Análisis Genómico (CNAG-CRG, Spain) using Illumina HiSeq. 2000 technology. A total of 56,327,604 paired-end, 126 bp sequence reads were generated. PrinseqQuality (http://prinseq.sourceforge.net/) was applied to quality filtering/trimming of reads (cut-off value, 20), and only reads with length ≥60-nt were used. Reads were assembled using the CLC Genomics Workbench version 5.0 (CLC Bio).
PacBio sequencing and de novo assembly. The single-molecule real-time (SMRT) sequencing technology developed by Pacific Biosciences (PacBio) 29 was used for long reads sequencing. A total of 311,471 pre-filtered reads were generated on a PacBio RS II sequencing instrument. The sequencing service was provided by the Norwegian Sequencing Centre (www.sequencing.uio.no), a national technology platform hosted by the University of Oslo and supported by the "Functional Genomics" and "Infrastructure" programs of the Research Council of Norway and the Southeastern Regional Health Authorities.
De novo genome assembly was carried out following a hierarchical genome-assembly process (HGAP), using the RS_HGAP_Assembly-3 protocol included in the PacBio SMRT Analysis Software v2.3.0 (Ref. 30   Assembly refinements. The 85 contigs, initially assembled by HGAP from the PacBio reads, were reanalyzed in order to discard those having a disproportionately low coverage (<40x) or short length (<15-Kb). Hence, 41 of those contigs were found to represent "spurious or artefactual" contigs, and were consequently discarded. After filtering, a total of 44 contigs were selected as bona fide genomic sequences. Twenty-eight of these contigs were found to correspond to complete chromosomes. For accurately assembling of the rest of chromosomes, different software packages and approaches were used 34 . Firstly, the Illumina paired sequencing reads were used to assess the possibility of joining some contigs by the SSPACE tool 35 . This approach allowed the complete assembly of chromosomes 7 and 35. The accuracy of these and the rest of assemblies was monitored by alignment of sequencing reads and its visualization by IGV 36 . Contigs belonging to the chromosomes 12, 15, 22, 26 and 28, were joined by the minimus2 assembler 37 , which uses an algorithm that calculates overlaps between contigs. On the other hand, the two contigs forming the chromosome 33 could be joined by means of the SSPACE-LongRead tool 38 , which selects and uses only the longest PacBio reads to construct a scaffolding. Finally, the gap size between pairs of contig was calculated (lower than 5-Kb in all cases) and closed with Gapfiller 39 , using the distance information derived from the paired-read data.
The contigs generated from the Illumina sequencing reads were aligned to the de novo assembled chromosomes using LAST aligner (http://last.cbrc.jp/). This allowed the identification of some Illumina contigs that aligned with the chromosomal ends but had overhanging sequences. In those cases, several tools were used to further extend the chromosomal ends. Thus, for chromosomes 4, 14 and 24, the optimal extension was attained with MAFFT multiple-aligner software 40 . For chromosomes 34 and 15, the best extension was obtained by minimus2. Finally, SSPACE-standard was found useful to extend a few nucleotides at the ends of the rest of chromosomes. A scheme of the complete pipeline is shown in Fig. 1.
Finally, a coverage analysis on the newly assembled chromosomes was performed using both Illumina and PacBio reads. Illumina reads were aligned by Bowtie2 (ref. 41 ), and PacBio bax.h5 reads were aligned by pbalign (which uses the BLASR method 42 ). Coverage analysis was done from each alignment along the 36 chromosomes using the GenomeCoverageBed tool (http://bedtools.readthedocs.io/en/latest/content/tools/genomecov.html). The graphical coverage plots files were generated with GNUPLOT (http://www.gnuplot.info/). Annotation of protein-coding genes and known non-coding RNAs. Bulk annotation of the assembled L. infantum genome was performed using Companion web server 43 , using the default settings and selecting the L. major (Friedlin strain) annotation as a reference. Given the importance of maintaining current L. infantum gene ID names, OrthoMCL 44 , BLAST searches were performed to assign correspondences between the Companion annotated genes and current L. infantum gene names (version 9, TriTrypDB.org). These additional data were combined with those provided by Companion annotation into a GFF3 file using an in-house Python script. Finally, the GFF3 file was manually curated to resolve ambiguous annotations and to name those new