Hybrid de novo whole-genome assembly and annotation of the model tapeworm Hymenolepis diminuta

Despite the use of Hymenolepis diminuta as a model organism in experimental parasitology, a full genome description has not yet been published. Here we present a hybrid de novo genome assembly based on complementary sequencing technologies and methods. The combination of Illumina paired-end, Illumina mate-pair and Oxford Nanopore Technology reads greatly improved the assembly of the H. diminuta genome. Our results indicate that the hybrid sequencing approach is the method of choice for obtaining high-quality data. The final genome assembly is 177 Mbp with contig N50 size of 75 kbp and a scaffold N50 size of 2.3 Mbp. We obtained one of the most complete cestode genome assemblies and annotated 15,169 potential protein-coding genes. The obtained data may help explain cestode gene function and better clarify the evolution of its gene families, and thus the adaptive features evolved during millennia of co-evolution with their hosts.


Background & Summary
The study of the genomics and transcriptomics of parasite model species has led to advances in the basic aspects of parasite biology, as well as new trends in human and veterinary medicine. Modern genomic tools, especially those based on a combination of multiple methods, allow detailed analyses of genome structure.
Our study used hybrid genome sequencing to examine the genome of the tapeworm Hymenolepis diminuta by three technologies: Illumina sequencing pair-end, Illumina mate-pair and MinION Oxford Nanopore DNA sequencing. H. diminuta is a well-described representative of the class Cestoda, the large group of parasitic flatworms that includes members known to be serious pathogens of vertebrate animals and humans 1,2 .
H. diminuta was chosen for the present study since it is commonly used in studies of new therapeutics, biochemical processes, immune responses and other host-parasite interrelationships during cestodiasis 3-8 and is considered the most important model species in experimental cestodology. Hymenolepis diminuta has a worldwide distribution as an intestinal parasite of rodents (primarily rats) and humans 9 , and the tapeworms of the genus Hymenolepis are considered to be among the most frequent causative agents of the human cestodiasis 1 .
Despite its importance as a model organism, the genome of H. diminuta is available only as a draft genome acquired as part of the 50 Helminth Genomes project initiative 10 . In addition, Gauci et al. 11 using the example of Echinococcus granulosus, highlight the possible limitations of published draft genomes of selected tapeworm species, one being the fact that they were sequenced using only Illumina short reads technology. Therefore, the ultimate goal of our study was to improve the accuracy of the draft genome of H. diminuta, by integrating data from www.nature.com/scientificdata www.nature.com/scientificdata/ three complementary approaches; this approach may significantly enhance the scientific value of the achieved datasets for future studies.
The combination of the recent progress in sequencing technologies and unlimited access to genomic data has fueled rapid development in the biomedical sciences, including parasitology. Most recently, the International Helminth Genomes Consortium released draft genomes (both published and unpublished) to present lineage-specific trends rather than individual species-specific differences 10 . This dataset of helminths genomes provides a number of new details important in studies of parasitic worms; however, there is an urgent need to continue helminth genome sequencing and improve the available genomes. An examination of the number and organization of EG95 E. granulosus vaccine-related encoding genes based on two available E. granulosus draft genomes published by Tsai et al. 12 and Zheng et al. 13 indicates that the genome sequence data available for E. granulosus offers limited potential for practical use 11 ; in part, it was not possible to map any of the EG95 gene family members previously characterized by cloning and sequencing genomic DNA fragments. These results have revealed limitations in available genomic data and highlights deficiencies present in current genomic resources, and thus, reinforced the need to supplement available datasets with new sequencing results. This can be achieved by simultaneous use of available sequencing technologies, providing both short and long reads. In recent years, such a hybrid approach has proven to be useful in improving quality of genome assemblies and improving discovery of gene family expansions. For instance, a hybrid approach was introduced for de novo human genome sequence 14 , one of the best described genomes. When assembling clownfish genome using high-coverage Illumina short reads and low-coverage Nanopore long reads, Tan et al. 15 observed substantial improvement in the genome statistics when compared with Illumina-only assembly. They suggest that development and improvement of Nanopore technology will shift toward the use of high-coverage long read-only assembly, followed by multiple iterations of genome polishing using Illumina reads. Genome improvements due to the use of hybrid sequencing have been applied to characterize the genetic polymorphism in Wuchereria bancrofti populations, and provide, among others, a list of genetic markers useful for monitoring changes in parasite genetic diversity 16 .
The present paper provides the first results of hybrid de novo whole-genome sequencing of H. diminuta combined with RNAseq analysis. Our assembly appears to be more complete than that available in WormBase ParaSite 17 and offers improved genome statistics. In this respect our results suggest that the procedure yielded one of the most comprehensive tapeworm genome assemblies available. In addition, our results are supported with RNA-seq analyses, which allow a better overview of the entire structure of the H. diminuta genome.
Here we confirm that the hybrid sequencing approach is the optimal method for obtaining the high quality data resulting in determination of a complete genome sequence. This cost-effective approach combining Illumina paired-end, mate-paired, and MinION Nanopore long reads allowed the retrieval of one of the most comprehensive tapeworm (or any parasitic worm) genome available, complimented by RNA sequencing data. These may result in better understanding of the biology of the parasite, its genetic diversity, adaptation to parasitic way of life and may allow new treatments and/or diagnostic tools to be identified in the near future. cultivation of H. diminuta adult cestodes. Six-week-old H. diminuta cysticercoids were removed from dissected Tenebrio molitor beetles under a microscope (100× magnification). Ten three-month-old rats were infected by voluntary oral uptake of six cysticercoids of H. diminuta per rat. Smears of their fecal samples were examined under a microscope (magnification 400×) five to six weeks from the initial infection, to verify the presence of adult parasites by their eggs. The rats were euthanized with 100 mg/kg intraperitoneal thiopental anaesthesia (Biochemie GmbH, Austria). The small intestines were removed immediately, adult parasites were isolated and washed up to 5× with 100 mM PBS with antibiotics added (1% penicillin) to remove debris. DNa isolation. Briefly after recovery from host intestine, DNA was isolated from tapeworm fragments containing only scolex and immature proglottids. Genomic DNA was isolated using a Genomic Midi AX isolation kit with ion-exchange membranes (A&A Biotechnology, Gdynia, Poland) according to the manufacturer's instructions. The integrity of the genomic DNA molecules was checked using agarose gel electrophoresis. The obtained DNA extracts were used immediately or stored at −20 °C until use.

Methods
RNa isolation and sequencing. A total of three adult H. diminuta tapeworms were homogenized in RLT buffer and total RNA was isolated from the homogenate using RNeasy Midi Kit (Qiagen, Germany). The sequencing library was prepared from 1 μg total RNA using TruSeq RNA Sample Preparation v2 Kit (Illumina, San Diego, CA, USA) according to manufacturer's instructions; the library was paired-end sequenced (2 × 100 bp) on the Illumina HiSeq 1500 platform.

WGS library preparation and sequencing.
For whole genome sequencing (WGS) 2.5 μg of high quality genomic DNA was used. Prior to the library preparation DNA was fragmented using Covaris M220 (Covaris, Inc, Woburn, MA, USA) and size selection was performed using BluePippin (Sage Science, Inc, Beverly, MA, USA) for the average insert size 600 bp. The library was prepared using NEBNext Ultra ® II DNA Library Prep Kit (New England BioLabs, Inc, Ipswich, MA, USA) according to manufacturer's instruction.
www.nature.com/scientificdata www.nature.com/scientificdata/ For mate-pair whole genome sequencing (MP-WGS) two different libraries, with (4 μg input DNA) and without (1 μg input DNA) size selection, were prepared. Libraries were constructed using Nextera Mate Pair Library Preparation Kit (Illumina) according to manufacturer's instruction. Size selection was performed using BluePippin (Sage Science) for fragments ranging from 5000 to 10000 bp (average size 8000 bp). The mean fragment size for the library without size selection was 2000 bp.
For Oxford Nanopore sequencing (ONT) high molecular DNA was isolated from tapeworm using phenol-chloroform extraction. Briefly, 200 mg of tapeworm tissue sample was washed twice with PBS buffer to remove excess rat stool material. After washing, the sample was submerged in 900 μl of TE buffer. The sample was lysed by the addition of 90 μl of 10% SDS, 10 μl of Proteinase K (20 mg/ml) and incubated at 37 °C for one hour www.nature.com/scientificdata www.nature.com/scientificdata/ until all cells were lysed. Following this, 200 μl of 5 M NaCl was added to the cleared lysate, which was subjected to phenol:chloroform:isoamyl alcohol extraction until no protein debris was visible in the interphase. After protein removal, the DNA was precipitated with isopropanol (0.7 volume added) and centrifuged for 10 minutes at 14000 rpm and washed with 70% ethanol. The DNA pellet was dried for a short time at room temperature and re-suspended in 100 μl of low-TE buffer (10 mM Tris and 0.1 mM EDTA pH = 8.0) containing RNase (50 μg/ml). DNA quality and integrity were checked using electrophoresis in standard 1% agarose gel and by PFGE using Biorad CHEF-II instrument. DNA quantity was measured with Qubit 3.0 fluorimeter and Broad Range chemistry (Thermo Scientific, Life Technologies).
The Oxford nanopore library was constructed by 1D ligation using two strategies. In the first, 8 μg of DNA was sheared into 20 kbp fragments using Covaris g-Tube and 5 μg of sheared template was taken for 1D library  www.nature.com/scientificdata www.nature.com/scientificdata/ construction using SQK-LSK108 kit (Oxford Nanopore Technologies). Approximately 1 μg of library was loaded into R9.4 flowcell system and sequenced on a MinION instrument for 24 hours. In the second approach, 20 μg of DNA was sheared into 20 kbp fragments followed by size selection on BluePippin instrument (Sage Science). Fragments above 10 kbp were recovered using PAC 30 kb cassette. 5 μg of recovered DNA was taken for 1D library construction using SQK-LSK108 kit and 1.5 μg of final library was loaded into R9.4.1 flowcell and sequenced on MinION sequencer.
De novo genome assembly. A hybrid assembly approach was employed, with several types of reads used in the assembly. Firstly, the datasets created from high quality of DNA reads from Illumina paired-end and Illumina mate-pair sequencing were assembled using tools based on de Bruijn graph, ABySS 18 and dnaasm 19 . The software versions are reported in Table 2. Secondly, the set of contigs (results of assembly) were combined based on the Oxford Nanopore long reads using two different tools: LINKS 20 and dnaasm-link 21 . This step was developed in an iterative way: firstly, results obtained from only short DNA reads were linked, where distance parameter in LINKS tool was set to 6 kbp. The obtained results were linked with those obtained for distance values of 7 kbp, then 8 kbp, 19 kbp, 20 kbp and 30 kbp. However, as the LINKS application requires a very large amount of RAM, the procedure   www.nature.com/scientificdata www.nature.com/scientificdata/ was performed using dnaasm-link running on C++ instead of Perl. In addition, dnaasm-link has a module to fill the gaps between contigs using sub-sequences from long DNA reads.
Functional annotation. The annotation pipeline was run using newly-obtained transcriptomic and genomic data from H. diminuta. During the first step, the RNA-seq data were mapped to the assembled genomic scaffolds using the STAR aligner. Obtained BAM file and genomic scaffolds were analyzed with BRAKER2 software with the Augustus tool to acquire the protein-expressing coding sequences. In the next step, BRAKER2 (amino acid sequences) and Trinity (transcriptomic sequences) outputs were used to obtain detailed genomic annotations using using MAKER2 pipeline (with -est2genome = 1; -prot2genome = 1).   www.nature.com/scientificdata www.nature.com/scientificdata/ All de novo assembled transcripts were searched against UniProt/SwissProt 22 database using BLASTx and BLASTp with an e-value < 10 −5 . Open reading frames (ORFs) were predicted using Transdecoder. The remaining functional annotation was obtained using g:Profiler and Trinotate pipeline, which uses several software packages: Hmmer, a protein domain identification (Pfam) tool, Rnammer to predicts ribosomal RNA and SignalP to predicts signal peptide sites.
Mitochondrial genome. Mitochondrial DNA was obtained and sequenced with Illumina technology as described above. The mtDNA was bioinformatically obtained from de novo assembly from the PET1 dataset (S59, S66, S70, S13, S41, S34 and S36 sets of reads). The mitogenome was analyzed and and characterized using CLC Main Workbench and MacVector software. The organization of mitochondrial genome is given in the 'Technical validation' section, where it is also compared with NC_002767.

Data Records
Data supporting the results of this article has been deposited at European Nucleotide Archive (EMBL). The study titled 'Hybrid sequencing of Hymenolepis diminuta genome' got Access Number ERP113437 23 , the project identifier is PRJEB30942. Raw Illumina and Nanopore reads have been given the indexes ERS3052629-ERS3052634, the assembly output is deposited under name 'H.diminuta_WMSil1' and identifier GCA_902177915 24 , mitochondrial genome under name 'Hymenolepis diminuta strain WMSil1 genome assembly, organelle: mitochondrion' , LR536429 25 . Annotation is included. Supporting data, also including script parameters, are available at figshare 26 .

technical Validation
Paired-end reads. Firstly, the quality of input data was checked using FastQC tool. The results confirmed the high quality of DNA reads -the reports were collected by the MultiQC tool and are available online at https:// doi.org/10.6084/m9.figshare.8798111.v1. Following this, the basic statistics of the paired-end tags were studied using the BBmap package (Table 1).
Further analysis used two data sets: PET1 and PET2. PET1 is a set with coverage 150× created from S59, S66, S70, S13, S41, S34 and S36, while PET2 was created from S3 with S47, and has 370× coverage. www.nature.com/scientificdata www.nature.com/scientificdata/ Mate-pairs reads. MP-WGS sequencing identified two sets of reads: MP1 and MP2. Both data sets consist of 100 bp DNA reads; however, MP1 was found to include 69, 558, 283 raw pairs of reads while MP2 had 54, 688, 723. Read quality, determined by FastQC, showed problems with adapter content, over-represented sequences and per sequence GC content, which is typical for this type of sequencing. To overcome this issue, the NxTrim 27 tool was used to filter only correctly paired reads based on the adapter location ( Table 3). The resultant sets of DNA reads were again checked by FastQC, and no such problems with DNA reads were observed.
After rejecting improperly paired DNA reads, the insert size value of the remaining mate pairs were examined using BBmap package. The results are presented in Fig. 1. ONt reads. Sequencing of nanopore library without size selection (ONT1) yielded 546,222 reads and 3.5 GB of sequence data with a mean read length of 6.37 kbp. Size selected library (ONT2) sequencing yielded 156,168 reads and 1.6 GB of sequence data with the mean read length of 10.1 kbp. Nanopore sequencing yielded totally 702,390 reads and 5.1 GB of data.
Raw nanopore data was base-called using Albacore (Oxford Nanopore Technologies, Oxford, UK). After quality filtering for quality and residual adapter removal using NanoFilt and Porechop. Long nanopore read data statistics, generated using NanoPlot, are presented in Figs. 2 and 3.
The error rate was checked using the BBmap package by mapping long DNA reads to contigs produced only from paired-end tags. The obtained results indicated a 25% error rate from nanopore DNA reads; therefore, the raw DNA reads were corrected using a Canu 28 correcting module, resulting in the error rate falling to 11%. Fig. 6 The organization of mitochondrial genome of Hymenolepis diminuta (WMS-il1 strain). All genes are transcribed in the same direction. The two leucine tRNA genes are designated by tRNA-LeuCUN and tRNA-LeuUUR, respectively, and two serine tRNA genes by tRNA-SerUCN and tRNA-SerAGN, respectively. Gene scaling is only approximate.
www.nature.com/scientificdata www.nature.com/scientificdata/ De novo assembly results. A hybrid assembly approach was employed, where short paired-end reads PET1 and PET2 datasets (depicted in Table 1) and mate-pair reads MP1 and MP2 (given in Table 3) and ONT1 and ONT2 long reads were used together.
The present study investigates the effect of applying reads from third-generation sequencers on de novo assembly results. In a typical de novo project, sequencing and assembly are performed iteratively until the results are of good enough quality and funds still remain. During each iteration, sufficient funds need to remain available for the next sequencing process, because the assembly costs are lower. From this point of view, two approaches can be used when performing a new experiment: (1) use the sequencing technology previously used in the project, or (2) complement results with sequencing technology not used previously in the project. Our results indicate that option (2) is a better choice, as adding results from new sequencing technology gives better statistics than additional reads obtained by the previously used technology.
As depicted in Table 4 we observed a significant improvement in assembly results between column 2 and column 3, when mate-pair reads were added, and between column 4 and column 5, when Nanopore reads were added. The improvements in assembly between column 2 and column 1, between column 4 and column 3, and  www.nature.com/scientificdata www.nature.com/scientificdata/ between column 6 and column 5 were less pronounced since they were obtained using the same sequencing methods. In particular, using N50 statistics as a measure of quality, we observed 21% better results (from 69.7 kbp to 84.2 kbp) when using PET1 + PET2 reads instead of PET1 reads (the sequencing coverage increases from 150× into 520×). Adding the mate-pair MP1 dataset (sequencing coverage 77×) into PET1 + PET2 (dataset has 597× coverage instead of 520×) improved N50 by 1000% (from 84.2 kbp to 842.2 kbp). The next mate-pair dataset, MP2 (sequencing coverage 61×, therefore all reads cover genome 658×) improved N50 by 0.2% (from 842.2 kbp to 844.2 kbp). Using the Nanopore dataset (ONT1, coverage 19×, mean read length 6.4 kbp) improved N50 by 206% (from 844.2 kbp to 1.7 Mbp), and the next Nanopore dataset ONT2 (coverage 8×, mean read length 10.1 kbp) improved N50 by 134%. A similar effect was observed when using a number of scaffolds.
In addition, our proposed approach is cost-and time-effective, and limited basically by the access to diverse sequencing technologies.
De novo transcriptome assembly. The Trimmomatic tool was used to trim out adaptors and low-quality fragments (Phred < 30) from the raw data. Reads shorter than 90 bp were removed from the dataset. Processed sequences were de novo assembled with Trinity with default parameters (k-mer = 25). This allowed to obtain a reference transcriptome comprising 28,282 transcripts. To confirm compatibility of RNA-Seq and DNA-Seq datasets, whole-transcriptiome mapping was performed to genomic scaffolds using BBMap, obtaining 85.65% (24,223/28,282) uniquely aligned transcripts.  Table 7. Overall statistics of the genome annotation. www.nature.com/scientificdata www.nature.com/scientificdata/ We used the BUSCO tool on the transciptome, yielding 784 complete, 668 complete and single-copy, 116 complete and duplicated, 40 fragmented and 154 missing BUSCOs. This result is better than the results of the scaffold analysis (Table 4).
Genome characteristics. Firstly, the k-mers distribution of the genome was studied using Jellyfish and GenomeScope tools. Jellyfish was used to obtain 51-mer count histogram in a subset of 7 GB of the short DNA reads, which was used to estimate genome size, heterozygosity and repeat content with the aid of GenomeScope. The size of the test genome was found to be approximately 185 Mbp (value close to the 177 Mbp resultant assembling size, see Table 5) with low heterozygosity (below 0.05%) and 15.4% repeat content (Fig. 4).
We tried to confirm the high repeat content value by launching the RepeatMasker 29,30 tool with Repbase 31 database (databases Dfam_Consensus-20170127 and RepBase-20181026). Several families of repeat elements covering only 0.72% of the genome were identified. However, in the presented genome assembly 9.127 Mbp of the 177.074 Mbp (5.2%) is known as 'N' signs. In addition, we estimate that approximately 8 Mbp (4.5%) of the genome has not been assembled. Most of the 'N' signs and unassembled sequences may consist of repetitive sequences, which may be a response to the high value of the predicted repeat content.
Genome assembly results comparison to results available at WormBase ParaSite. The H. diminuta genome has previously been studied and the genome draft is available 10 . However, our sequencing effort 23 resulted in approximately 45× better N50 statistics (2.3 Mbp versus 51.0 kbp; in presented study 13× fewer scaffolds were obtained: 719 scaffolds in comparison with 9867, with the longest being almost 7 Mbp compared to 356 kbp in the previous work (Table 5). Our results were also evaluated using the Circoletto tool; example results are presented in Fig. 5.
The rrnL gene (967 bp) is separated from the rrnS gene (709 bp) by the tRNA-Cys gene. The length of the tRNA genes vary from 59 bp (tRNA-Ser) to 72 bp (tRNA-His). The 12 protein-coding genes encoded a total number of 3,363 amino acids. The total length of all protein-coding genes was found to be 10,089 bp. The length of the individual protein-coding genes varied from 261 bp (nad4L gene) to 1599 bp (cox1 gene). Except for the nad4 gene, all the protein-coding genes use the ATG start codon, whereas the nad4 gene uses ATT as a start codon. The majority www.nature.com/scientificdata www.nature.com/scientificdata/ of identified protein-coding genes are terminated with the TAG termination codon; the only exceptions are cox2 and nad6 genes, which are terminated with the TAA codon.
In the mitogenome of H. diminuta (WMS-il1 strain) two non-coding regions were found: the larger between ND5 and tRNA-Gly genes, and the shorter is between the tRNA-Tyr and tRNA-Ser genes. The nucleotide composition of the obtained mitogenome is A = 25.4%, T = 45.6%, G = 19.3% and C = 9.6%.
Our mitogenome analysis of the H. diminuta WMS-il1 strain mitogenome was performed using data from the Illumina next-generation sequencing. All 36 genes previously found in mitogenomes of other cestode species were identified [32][33][34] . The length, structure and composition of the coding regions are also similar to these previously described in tapeworms, including NCBI H. diminuta reference sequence NC_002767.1 35 . No differences were observed in the gene sequence encoding tRNA. However, both rRNA-coding genes differed with regard to two bases when compared to the reference sequence. Interestingly, the protein-coding regions showed substantial variability and only ND3 was identical as these described in reference sequence. These differences are shown in Table 6.
Gene prediction. Gene prediction was performed with genomic scaffolds according to the protein sequences of H. diminuta (PRJEB507) and other closely-related organisms: H. nana (PRJEB508), H. microstoma (PRJEB124) and Echinococcus multilocularis (PRJEB122), downloaded from WormBase ParaSite database 36 Version: WBPS12 (WS267). This step was processed again by MAKER2 software (with -est2genome = 0; -prot2genome = 0). The annotation files (GFF3) obtained from each species were combined and both results were compared using custom script in the R environment ver. 3.5.0. Next, CDS annotations not confirmed in either pathway which were shorter than 150 nt (as suggested by NCBI) were removed from the final GFF3 file using Genome Annotation Generator (GAG) 37 with -rcs 150 option. The general statistics of GFF file modifications using GAG are presented in Table 7.
A total of 15,169 potential protein-coding genes were predicted in the assembled H. diminuta genome and functionally annotated, which encodes 19,651 mRNAs. For extracting CDS sequences, the gffread (https://github. com/gpertea/gffread) script was applied. In total 16,983 (86.42%) homologs were identified in H. diminuta with a median sequence identity of 98.91%, 15,144 (77.06%) homologs in H. microstoma with a median sequence identity of 80.36%, 14,668 (74.74%) homologs in H. nana with a median sequence identity of 78.04%, and 14,132 (71.91%) homologs in E. multilocularis with a median sequence identity of 60.00% ( Fig. 7(a,b)) by searching WormBase ParaSite database using BLASTp 38 and CDS sequences as query.
annotation results. Our sequencing and annotation results enriched de novo assembly reference of the H.diminuta genome available from WormBase ParaSite. A considerable body of the annotation created only by in silico prediction is incomplete, and requires re-annotation. Our acquisition of RNA-seq data offers a significant improvement in the finalization of the annotation processes, as even mRNA sequences from related organisms (H. microstoma, H. nana, E. multilocularis) do not always form the best basis for exon-intron structure prediction. By using transcriptome evidence from the same species (H. diminuta) it was possible to confirm intronic donor-acceptor sites according to the alignment of cDNA and genomic DNA. Our improved annotation allowed the splice site to be corrected according to de novo assembly transcriptome aligned to H. diminuta genome. Software applied in this study allowed us to add UTR regions to previously-annotated genes ( Fig. 8a and Suppl. A). Our data includes some fixes of the reference CDS regions (Fig. 8b and Suppl. B); in addition, the H. diminuta genome was supplemented with genes that have not yet been annotated in the reference genome ( Fig. 8c and Suppl. C). In some cases, two gene annotations, predicted by Sanger Institute (annotated on two separated scaffolds, blue -HDID_scaffold0000291 and orange -HDID_scaffold0000029 bars on the Fig. 8d and Suppl. D) were merged into individual complete protein-coding gene.

code availability
The software packages used for the analysis with version numbers, are given in Table 2. The parameters of scripts used for the analysis are available at figshare 26 . www.nature.com/scientificdata www.nature.com/scientificdata/ www.nature.com/scientificdata www.nature.com/scientificdata/ author contributions