A chromosome-scale reference genome for Giardia intestinalis WB

Giardia intestinalis is a protist causing diarrhea in humans. The first G. intestinalis genome, from the WB isolate, was published more than ten years ago, and has been widely used as the reference genome for Giardia research. However, the genome is fragmented, thus hindering research at the chromosomal level. We re-sequenced the Giardia genome with Pacbio long-read sequencing technology and obtained a new reference genome, which was assembled into near-complete chromosomes with only four internal gaps at long repeats. This new genome is not only more complete but also better annotated at both structural and functional levels, providing more details about gene families, gene organizations and chromosomal structure. This near-complete reference genome will be a valuable resource for the Giardia community and protist research. It also showcases how a fragmented genome can be improved with long-read sequencing technology completed with optical maps.


Background & Summary
Giardia intestinalis has a duplicated cell structure with two transcriptionally active nuclei. It is one of the most prevalent parasitic protists causing diarrhea, infecting a wide range of hosts, including humans. Giardia has been categorized into eight different assemblages (genotypes, A-H) based on host-specificity and genetic differences 1 . Assemblage A and B infect humans, with WB isolate from assemblage A being the most extensively studied. The first sequenced G. intestinalis genome is from the WB isolate, and was published in 2007 2 . At the time, 200,000 reads from both ends of small-insert plasmid libraries and 2,400 end sequences from a 200 kbp library in bacterial artificial chromosome (BAC) vectors were generated and sequenced with Sanger sequencing. This yielded a genome of 11.7 Mbp in size, distributed on 306 contigs (92 scaffolds) 2 (Table 1). Since then, this genome has been used as the reference genome for Giardia research. Subsequently, a few other Giardia genomes have been published, including human isolates (GS, GS-B and BAH15c1) from assemblage B [3][4][5] , human isolates (DH, AS175 and AS98) from assemblage A2 4,6 , one pig isolate (P15) from assemblage E 7 and two dog isolates from assemblages C and D 8 . The availability of several genomes has advanced Giardia research and it is now possible to use different omics techniques during studies in vitro and in vivo.
As part of diplomonads, a group of unicellular protists, Giardia genomes have made it possible to perform comparative genomic studies within diplomonads. Comparative genomics between the genomes of the salmon parasite Spironucleus salmonicida and Giardia revealed how S. salmonicida adapted to colonize different sites in the host 9 . Comparative studies between the transcriptome of free-living diplomonad Trepomonas and parasitic diplomonad genomes revealed how a free-living organism evolved from its parasitic ancestor 10 .
However, all published Giardia genomes are fragmented, due to the limitation of the sequencing technologies used at the time. Fragmented genomes have limited our understanding on chromosome structure and how the genome evolves at the chromosomal level. During the last years, long-read third generation sequencing has matured. For example, Pacbio claims that its newest HiFi sequencing can generate reads up to 300 kbp with half of the bases in reads >160 kbp 11 .
To obtain a high-quality reference genome of G. intestinalis, we have re-sequenced the WB genome using Pacbio technology. Long-reads were assembled into long contigs, which were then scaffolded into near-complete five chromosomes with the aid of optical maps 12 . Using a comparative genomics approach and RNA-Seq reads mapping, we could also improve both the structural and functional annotations.

Methods
DNa preparation and sequencing. G. intestinalis WB/C6 (ATCC 50803) was cultivated according to ref. 13 in TYDK media in tightly capped slanted culture tubes at 37 °C. Total genomic DNA was extracted from 2.3 × 10 8 trophozoites using Qiagen Blood & Tissue Kit 50 (per manufacturer's instructions), and purified using phenol-chloroform extraction. It was then further purified with Qiagen Genomic-tip 20/G. The concentration and quality of the extracted DNA was determined by NanoDrop and agrose gel electrophoresis. Genomic DNA (5.6 μg) was sequenced with PacBio RS II using 8 SMRT cells at the Uppsala Genome Center at the Science for Life Laboratory (Uppsala University), which generated 411,835 reads in 3.6 billion bases (Table 1). Reads have an N50 length of 12.7 kbp and 10% them are longer than 20 kbp.
To assist base error correction, the PacBIO long-read data was supplemented with short-read (Illumina) sequence data at a high coverage (100X). These short-read data was generated from G. intestinalis WB/C6 grown in TYDK media as per 13 in a 25 mL tissue culture flask. Total genomic DNA was extracted from 1 × 10 7 trophozoites using the Qiagen Powersoil Kit (per the manufacturer's instructions). Approximately 1 μg of genomic DNA was prepared for Illumina Truseq library construction (fragment size 300 bp) and subjected to paired-end (75 bp) sequencing on Illumina HiSeq 2000 per the manufacturer's instructions. This generated 16 million paired-end reads.
Software mentioned below without version information were all from the SMRT Analysis pipeline. Reads were assembled de novo with HGAP 14 followed by consensus sequence calling with Quiver 14 . This yielded 80 contigs. Those contigs were then mapped to the optical maps of the five chromosomes (Fig. 1a) 12 using MapSolver (v3.2.0) provided by OpGen. Neighboring contigs were then stitched together as follows: If overlaps were identified between the neighboring contigs, the two contigs were merged based on the aligned overlapping sequence; If no overlap was detected, the two contigs were stitched with Ns in between to represent gap with size determined by the aligned optical map. The resulting five scaffolds were fed into BridgeMapper workflow 14 , which produces split alignments using BLASR. The split alignments were visualized using SMRT View, and no misassemblies and structural variation were detected. PBJelly (v15.8.24) 15 was used to further close the gaps in the five scaffolds using the Pacbio reads. The scaffolds were then further polished with Quiver 14 , and the reads failed to map to the scaffolds were assembled independently with canu (v1.4) 16 . Canu contigs and the five scaffolds were combined, and further polished with Quiver. This resulted in 35 scaffolds, with the five major ones representing 97% of the total size at  Table 1. Comparison of the old and the new G. intestinalis WB genome. *200,000 from end sequences with small insert, 2,400 end sequences from 200 kbp long insert. **There were 4 identified intron-containing genes in the first published draft genome. Not all the discovered intron-containing genes were properly integrated in the newest GiardiaDB either. This version integrates all the discovered introns, which is consistent with our results searching for the de novo introns.
12.6 Mbp, slightly bigger than the old draft genome ( Table 1). The five major scaffolds representing the five Giardia chromosomes contain only four internal gaps compared to the 137 gaps in the old version (Table 1, Fig. 1b). During the annotation stage below, we noticed that some full length genes were truncated in the new draft genome due to frame-shifts caused by homopolymer errors introduced by PacBio reads. Illumina DNA sequencing reads were mapped to the draft reference using BWA (v0.7.15) short read alignment component 17 and the alignment bam file was sorted by samtools (v1.8) 18 and its mpileup function was used to generate pileup file listing per base mapping information. Pileup file was parsed and indels and SNPs were updated when the alternative base/indel has at least 50% support from at least 10 reads mapped. This resulted in 50 SNPs updated affecting 20 genes, and 85 indels updated affecting 38 genes.
Heterozygosity estimation. Illumina reads were re-aligned to the base-corrected reference genome using BWA (v0.7.15) and samtools (v1.8) mpileup result was re-generated. SNP (or allelic sequence heterozygosity (ASH)) sites mentioned below were called in positions of at least 20X base coverage with an alternative base in at least 10% of the reads. Genome annotation. The previously published G. intestinalis WB genome was downloaded from GiardiaDB (v36.0) 19 and the annotation was transferred to the new reference genome using RATT (v0.95) 20 . GlimmerHMM (v3.0.1) 21 and Prodigal (v2.6.3) 22 were used for de novo gene prediction, and 500 RATT transferred annotations were selected to train GlimmerHMM prediction. Similarity information from BLASTP 23 against NR database and domain information from Conserved Domain (CD) search 24 were combined for functional annotation. RATT transferred annotations and de novo predictions were merged, and the ones inconsistent of the three sources were manually examined. RNA-Seq reads (SRR10063826 25 ) were mapped to the genome with BWA (v0.7.15), and the mapping information were used as a guideline for structural annotation. Manually curated annotations provided by Giardia researchers 19 were also incorporated into the functional annotation. Inconsistent annotations between the old version and the new version were double checked to ensure the updated annotation was an improvement. For multi-copy genes, synteny between the two versions were examined to determine the most appropriate geneid for the gene, an attempt to keep the synteny consistent between the two versions of genomes.  www.nature.com/scientificdata www.nature.com/scientificdata/ Introns in Giardia have been published in different papers 2,26,27 , but were not curated and integrated in the genome. In this version of annotation, we integrated all the 8 cis-spliced introns and 5 trans-spliced introns ( Table 1). Searching for the new general intron motif did not reveal new intron candidates.
We have in total annotated 4,963 protein-coding genes plus an additional 306 pseudo genes in the new genome (Table 1). Although the number of genes is reduced compared to the old genome, the new genome has on average longer genes and smaller intergenic regions ( Table 1). The coding percentage remains the same.
18S, 28S and 5S ribosomal RNAs (rRNAs) were predicted using rnammer (v1.2) 28 . 5.8S rRNA sequence was retrieved from NCBI, then searched against the new genome with BLASTN to locate all the copies. tRNAs were identified using tRNAscan-SE (v1.23) 29 . There are more duplicated copies of tRNAs and rRNAs in the new genome (Table 1). 5 S rRNAs are found internally, while 18S, 5.8S and 28S rRNAs are found at most of the chromosome ends (Fig. 2). Two partial 18S rRNAs sit at the two ends of the biggest internal gap at the 5′ end of chromosome 5 (Fig. 2), which suggest the unfilled gap could be caused by tandem copies of rRNAs. If so, it will make this ribosomal RNAs (28S + 5.8S + 18S) region 318 kbp long, approximately 86 copies of rRNAs.

Gene family annotation and distribution.
Giardia has a large group of cysteine-rich proteins, which are divided into 3 subgroups, including variant surface proteins (VSPs), high cysteine membrane proteins (HCMPs) and high cysteine proteins (HCPs). VSPs are secreted to the membrane and thus carry 5′ signal peptide and have a defined 3′ tail with a transmembrane domain followed by CRGKA pentapeptide motif. VSPs without proper signal peptide predicted by signalP (v5.0b) 30 were either adjusted to the right start codon, if possible, or annotated as pseudoVSPs (pVSP, ψVSP). On top of lacking proper signal peptide, ψVSPs also tend to sit in arrays without proper start codons. In the end, we have 133 VSPs compared to 196 in the old genome, and an additional 208 ψVSPs.
Among the 133 VSPs, there are 38 in duplicated pairs. 12 pairs have a → ← orientation, and 7 pairs have a ← → orientation. Among the 7 pairs pointing away from each other, 4 pairs contain a reverse transcriptase remnant in between, and one contains a NEK kinase in between ( Table 2). It was already observed, while comparing genomes from different Giardia isolates, that conserved synteny in Giardia genomes breaks at regions enriched with multi-gene families including cysteine-rich protein families, NEKs and ARPs 3 . With this near-complete genome, we could see how these gene families are distributed along the chromosomes (Fig. 2). The gene families are all over the chromosomes. Giardia genome is packed with genes with 81% genome coding and small intergenic region. However, regions around VSPs tends to be gene-poor (Fig. 2, green dots in coding% track), and those regions also have higher allelic sequence variation in the mapped reads (Fig. 2, red dots in SNPs track).

Data Records
This Whole Genome Shotgun project has been deposited at DDBJ/ENA/GenBank under the accession AACB00000000. The version described in this paper is version AACB03000000 33 , and the GenBank assembly accession is GCA_000002435.2 34 . To distinguish the duplicated gene from the original one, 'd' was added after locus_tag GL50803_; To distinguish this version of annotation from the old one, '00' was added after locus_ tag GL50803_ in all the new geneids. The datasets generated and analyzed in this project are also available via GiardiaDB. Raw DNA sequence reads from Pacbio and Illumina are deposited at NCBI Sequence Read Archive (SRA) under accession number SRP191500 35 . technical Validation Genome completeness. The five assembled chromosomes correlate well with the optical maps (Fig. 1a), with 92.2% of the optical maps covered by the assembled sequence. The biggest missing pieces lie at the 5′ end of the chromosome 5 (Fig. 1). From the sequences mapped to that region, we could see approximately four copies of non-LTR LINE-like retrotransposon GilM (Fig. 2) 36 . One copy of GilM is 5,474 bp long and encodes a reverse transcriptase (RT). Mapping all the raw reads to the region gives an average coverage at 1,800X. While the whole genome was sequenced at 200X coverage, it indicates that there are at least 9 copies of this GilM in the genome. The size and the repetitive pattern of the optical map at the 5′ end of the chromosome 5 (Fig. 1a) indicate more than 9 copies of GilM. This difference could be caused by subtelomeric size variation, which has been widely observed in humans and yeast 37,38 , and it has been observed that the chromosome sizes vary in A isolates 6 . Besides that, the chromosomal ends could have been sequenced at less coverage than average. Synteny between the old and the new genome is highly conserved, with the new genome being more complete. Four out of ten chromosome ends were assembled into telomeric repeats (TAGGG)n (Fig. 1b). The orientation of the chromosomal sequences in the new genome is ordered after the orientation of the optical maps, which is reverse complement to the old one (Fig. 1b).  www.nature.com/scientificdata www.nature.com/scientificdata/ The remaining 33 contigs are all small with sizes <33 kbp, and enriched with repeats, ribosomal RNAs and gene families.
The new genome is also more complete in reconstructing duplicated regions. In fact 6.9% of the genome resulted from duplications. Duplications are more commonly observed in short stretches of DNAs with just one gene, but multi-gene large duplications up to 20 kbp are also observed. This might be the largest duplication the assembly could resolve due to the limitation of read length. Members of multi-gene families are often seen in the duplicated regions. annotation improvement. The first published Giardia WB genome was over-annotated with small open reading frames (ORFs) overlapping real genes, as well as genes with wrongly assigned start codons (Fig. 3a). This is problematic for many analyses like unique genes, intergenic regions and regulatory elements. We used information from domain searches, start codon motifs, RNA-Seq read alignments and homologs from other Giardia genomes, to clean up and improve the structural and functional annotations.
Among the 4,963 protein-coding genes, 2,948 genes stay exactly the same, 17 stay the same but with SNPs, 1,173 genes have the same sequences but updated descriptions. 626 genes have adjusted start codons, based on the alignments of orthologous genes within Giardia as well as the motifs around the start codons. In fact, the new start codons often sit after a stretch of A-rich sequences (Fig. 3b), while this motif is lacking at the start codons of the old genes (Fig. 3c). This A-rich motif is also inline with the general sequence profile at the start codons for the whole genome (Fig. 3d) 2 . 55 genes are new to the new genome, and most of their sequences were found in the old genome but not annotated due to misannotations. 138 genes are new to the new genome but share sequence similarity with genes already annotated. 6 genes are the same but have part of the sequences changed due to frame-shifts in the new reference genome. Figure 3 showcases an 8 kbp region on chromosome 5 where the new genome has better annotation. In this 8-gene region, 3 of the genes were shortened with now the start codons better matching the RNA-Seq reads mapped (Fig. 3a), as well as the genes start with A-rich motif right next to or close to the start codons (Fig. 3b). All these three genes got updated functional annotations including domain information. Shortening GL50803_16996 also allowed space for a new gene GL50803_0060288, BolA-like protein. The two grey colored hypothetical genes in the old genome, overlapping completely with the functional genes, are most likely misannotations and were thus removed in the new genome.
The new genome has also improved annotation on hypothetical proteins with now 2,099 hypothetical proteins instead of 3,545 in the old genome. 955 hypothetical proteins in the old genome were removed in the new genome due to reasons stated above, and 69 hypothetical proteins are new to the new genome. 598 hypothetical proteins have new functional annotation based on domain and BLAST matches, and 38 functionally annotated genes in the old genome were changed to hypothetical proteins due to lack of support for the functions.