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 differences1. 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 20072. 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 B3,4,5, human isolates (DH, AS175 and AS98) from assemblage A24,6, one pig isolate (P15) from assemblage E7 and two dog isolates from assemblages C and D8. 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.

Table 1 Comparison of the old and the new G. intestinalis WB genome.
Full size table

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 host9. Comparative studies between the transcriptome of free-living diplomonad Trepomonas and parasitic diplomonad genomes revealed how a free-living organism evolved from its parasitic ancestor10.

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 kbp11.

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 maps12. 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 × 108 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 per13 in a 25 mL tissue culture flask. Total genomic DNA was extracted from 1 × 107 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.

Genome assembly

SMRT Analysis (v2.3.0) pipeline14 provided by Pacbio was used for the Pacbio reads. Software mentioned below without version information were all from the SMRT Analysis pipeline. Reads were assembled de novo with HGAP14 followed by consensus sequence calling with Quiver14. 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 workflow14, 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 Quiver14, 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 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).

Fig. 1
figure 1

Near-complete five chromosomes. (a) Restriction enzyme (MluI) maps of the five chromosomes aligned with the genomic sequences digested with MluI in silico. Each vertical line inside boxes represents a restriction enzyme cutting site. Gaps in the genomic sequences are represented with a horizontal line outside of boxes. (b) Circular plot comparing the old five chromosomes (to left) to the new five chromosomes (to right). Chromosomal sequences are represented in grey at the outermost circle with gaps in white bands and telomeric repeats in red. BLASTN matches between the two genomes are shown as blue ribbons in the middle. R package circlize (v0.4.8) was used to draw the circular plot39.

Full size image

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 component17 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 BLASTP23 against NR database and domain information from Conserved Domain (CD) search24 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 (SRR1006382625) 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 researchers19 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.

Introns in Giardia have been published in different papers2,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.

Fig. 2
figure 2

Circular plot of the five chromosomes. Chromosomal sequences are represented in grey at the outermost circle with gaps in white bands and telomeres in red. Inner tracks are arranged as: GC%, ARPs/NEKs, (ψ)Reverse transcriptases/rRNAs, (ψ)VSPs/HCMPs, Coding density, SNPs density, regions with similarity. Regions with similarity represent BLASTN matches against itself with >95% sequence identity and >2000 bp in length. The circular plot was drawn with R package circlize (v0.4.8)39.

Full size image

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).

Table 2 Arrangement of VSP pairs.
Full size table

NEK kinase and the arbitrarily named protein 21.1 share ankyrin-repeat domains, and there are 179 NEK kinases and 242 protein 21.1 in the old genome. In the new genome, ankyrin-repeat containing genes are named as follows: Proteins with both kinase domain and ankyrin-repeat are called NEK kinases, and there are 184 of them; Proteins with only ankyrin repeats are called ankyrin repeat protein 1 (ARP-1), and there are 267 (+5 pseudo copies) of them; Proteins with ankyrin repeat as well as zinc-finger domain(s) are called ankyrin repeat protein 2 (ARP-2), and there are 33 of them; Proteins with ankyrin repeat as well as domain(s) other than a zinc-finger domain are called ankyrin repeat protein 3 (ARP-3), and there are 5 of them. The domains were found using HMMER3 (v3.0)31 search against Pfam (v31.0)32.

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 ARPs3. 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 AACB0300000033, and the GenBank assembly accession is GCA_000002435.234. 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 SRP19150035.

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 yeast37,38, and it has been observed that the chromosome sizes vary in A isolates6. 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).

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.

Fig. 3
figure 3

An 8 kbp genomic region on chromosome 5 with improved annotation and start codon sequence logos. (a) An 8 kbp long genomic region located on chromosome 5. The first part shows the RNA-Seq reads mapped to the region. A coverage cutoff 30 was used for better display. The new genome is drawn directly below the RNA-Seq reads with the old genome aligned at the bottom. Orthologs between the two genomes are linked with light pink bands. Genes are shown in arrowed boxes filled with different colors. White indicates genes without modification, purple (60288) indicates the new BolA-like gene to the new genome, grey represents the unique genes to the old genome (deleted from the new genome). Blue represents genes with adjusted start codons, and in these three genes, they were shortened with updated descriptions. (b) Sequence logo at the start codon of the 626 new genes with updated start codons. (c) Sequence logo at the start codon of the 626 old genes from the old genome. (d) Sequence logo at the start codon of all the 4,963 protein-coding genes. Sequence logos were drawn with R package seqLogo (v1.50)40.

Full size image

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.