No evidence for extensive horizontal gene transfer in the genome of the tardigrade Hypsibius dujardini

Assembly of the Genome of H. dujardini.

Using propidium iodide flow cytometry, we estimated the genome of H. dujardini to be ∼110 Mb, similar to a previous estimate (20). We sequenced and assembled the genome of H. dujardini using Illumina short-read technology. Detailed methods are given in SI Appendix. Despite careful cleaning before extraction, genomic DNA samples of H. dujardini were contaminated with other taxa. Adult H. dujardini have only ∼10³ cells, and thus a very small mass of bacteria would yield equivalent representation in raw sequence data. A preliminary assembly (called nHd.1.0) generated for the purpose of contamination estimation spanned 185.8 Mb. We expected assembly components deriving from the H. dujardini genome to have similar proportion of G + C bases (GC%), and to have the same coverage in the raw data (because each segment is represented equally in every cell of the organism). Contaminants may have different average GC% and need not have the same coverage as true nuclear genome components. Taxon-annotated GC-coverage plots (TAGC plots or blobplots) (38, 39) were used to visualize the genome assembly and permitted identification of at least five distinct blobs of likely contaminant data with GC% and coverage distinct from the majority tardigrade sequence (Fig. 1A). Read pairs contributing to contigs with bacterial identification and no mitigating evidence of tardigrade-like properties (GC%, read coverage, and association with eukaryote-like sequences) were conservatively removed. There was minimal contamination with C. reinhardtii, the food source (41). Further rounds of assembly and blobplot analyses identified additional contaminant data (39), which was also removed. An optimized assembly, nHd.2.3, was made from the cleaned read set. Contigs and scaffolds below 500 bp were removed. Mapping of H. dujardini poly(A)+ mRNA-Seq (42) and transcriptome (12) data were equivalent between nHd.1.0 and nHd.2.3 (SI Appendix, Table S1); therefore, we conclude that we had not overcleaned the assembly.

• Download figure
• Open in new tab
• Download powerpoint

Fig. 1.

H. dujardini genome assembly. (A) Blobplot of the initial nHd.1.0 assembly, identifying significant contamination with a variety of bacterial genomes. Each scaffold is plotted based on its GC content (x axis) and coverage (y axis), with a diameter proportional to its length and colored by its assignment to phylum. The histograms above and to the right of the main plot sum contig spans for GC proportion bins and coverage bins, respectively. (B) Blobplot of the nHd.2.3 assembly (as in A). (C) Blobplot of the nHd.2.3 assembly, with scaffold points plotted as in B but colored by average base coverage from mapping of RNA-Seq data (42). A high-resolution version of this figure is available in SI Appendix.

The nHd.2.3 assembly had a span of 135 Mb, with an N50 length of 50.5 kb (Table 1). The assembly was judged relatively complete. It had good representation of a set of highly conserved, single-copy eukaryotic genes from the Core Eukaryotic Genes Mapping Approach (CEGMA) set (43), and these had a low duplication rate (1.3–1.5). A high proportion of H. dujardini mRNA-Seq (Fig. 1C) (42), transcriptome assembly (12), expressed sequence tags (ESTs), and genome survey sequences (GSSs) mapped to the assembly. We predicted a high-confidence set of 23,021 protein-coding genes using AUGUSTUS (44). The number of genes may be inflated because of fragmentation of the assembly, as 2,651 proteins lacked an initiation methionine, likely because they were at the ends of scaffolds, and were themselves short.

Assembly of the H. dujardini genome was not a simple task, and the nHd.2.3 assembly is likely to still contain contamination. We identified 327 scaffolds (5.0 Mb) that had read coverage similar to bona fide tardigrade scaffolds but similarity matches to bacterial genomes. Some of these scaffolds also encoded eukaryote-like genes and may represent HGT or misassemblies. Some scaffolds (195 spanning 1.5 Mb) had only bacterial or no genes and were very likely to be contamination. We identified no scaffolds with matches to bacterial ribosomal RNAs (rRNAs) but did find an 11-kb scaffold with best matches to rRNAs from bodonid kinetoplastid protozoa. Two additional small scaffolds (6 and 1 kb) encoded kinetoplastid genes (a retrotransposon and histone H2A, respectively). No other genes were found on these scaffolds, and their high coverage likely resulted from the loci being multicopy in the source genome. The genome was made openly available to browse and download on a BADGER (45) server at www.tardigrades.org in April 2014 (SI Appendix, Fig. S3).

Claims of Extensive Functional Horizontal Gene Transfer into H. dujardini.

Boothby et al. (13) published an estimate of the genome of H. dujardini, referred to as the UNC (University of North Carolina) assembly hereafter, based on a subculture of the same culture sampled for nHd.2.3. They suggested that the H. dujardini genome was 252 Mb in span, that the tardigrade had 39,532 protein coding genes, and that more than 17% of these genes (6,663) had been derived from extensive fHGT from a range of prokaryotic and microbial eukaryotic sources. Given this claim, and the striking difference between the UNC assembly and our assembly, we set out to test the hypothesis that these “HGTs” were in fact unrecognized contamination in the UNC assembly.

Surprisingly, the UNC assembly had poorer metrics than nHd.2.3 (31) (Table 1), despite the application of two independent long read technologies [Pacific Biosciences (PacBio) and Moleculo] and equivalent short read data. Scaffold N50 length was one-third that of nHd.2.3, despite UNC having discarded all scaffolds shorter than 2 kb. The UNC assembly span was 1.9 times that of nHd.2.3, in conflict with the UNC authors’ own (20) and our genome size estimates. The UNC protein prediction set was 1.7 times as large as that from nHd.2.3. UNC had good representation of CEGMA genes (Table 1), but contained more than three copies on average of each single-copy locus. Such multiplicity of representation of CEGMA single-copy genes can arise through bacterial contamination (as the CEGMA gene set is not explicitly designed to exclude loci with bacterial homologs; SI Appendix).

About one-third of the span of UNC does not appear to be derived from the tardigrade. Many scaffolds had low coverage compared with bona fide tardigrade scaffolds (Fig. 2A), had different relative coverages in different libraries (SI Appendix, Fig. S4 B–D), were not represented in our raw data (Fig. 2B), and had overwhelmingly noneukaryote taxonomic assignments (SI Appendix). The absence of all but marginal similarity to metazoan sequence also suggests that these contigs are not chimeric coassemblies. All of the longest scaffolds in UNC were bacterial (Fig. 3A), and few bacterial scaffolds had read coverage support in both UNC and Edinburgh raw data (Fig. 3B). We identified 15 scaffolds in UNC with high-identity matches to rRNA genes from Armatimonadetes, Bacteroidetes, Chloroflexi, Planctomycetes, Proteobacteria, and Verrucomicrobia (SI Appendix, Table S2). We also identified contamination that is likely to derive from other genomes. Two very similar UNC scaffolds (scaffold2445 and scaffold2691) both contained two tandemly repeated copies of the rRNAs of a bdelloid rotifer related to Adineta vaga. We found a large number of additional matches to the A. vaga genome (28) in UNC, but these may be bacterial contaminants matching A. vaga bacterially derived fHGT genes (28, 30). A total of 0.5 Mb of scaffolds had best sum matches to Rotifera rather than to any bacterial source (SI Appendix). Six mimiviral-like proteins were identified, five of which involved homologs of the same protein family (with domain of unknown function DUF2828). Mimiviruses are well known for their acquisition of foreign genes (46), and these scaffolds may derive from a mimivirus rather than the tardigrade genome. Overall, very few of the total fHGT candidates proposed by Boothby et al. (13) were in scaffolds that were not obviously contaminant (SI Appendix).

• Download figure
• Open in new tab
• Download powerpoint

Fig. 2.

Contaminants in the UNC assembly. (A) Blobplot of the UNC assembly with coverage derived from pooled UNC raw genomic data. (B) Blobplot showing the UNC assembly with coverage derived from the Edinburgh short insert genomic data. (C) Blobplot (as in A) with the scaffold points colored by average RNA-Seq base coverage. A high-resolution version of this figure is available in SI Appendix.

• Download figure
• Open in new tab
• Download powerpoint

Fig. 3.

Identifying HGT candidates. (A) Stacked histogram showing scaffolds assigned to different kingdoms (Bacteria, Eukaryota, and “no hits”) in different length classes for UNC and nHd.2.3 assemblies. The nHd.2.3 assembly had no scaffolds >1 Mb, and all of the longest scaffolds (>0.5 Mb) in the UNC assembly were bacterial. (B) Coverage-coverage plot of the UNC assembly using the Edinburgh short insert data (x axis) and in the pooled UNC short insert data (y axis). (C) Coverage-coverage plot of the nHd.2.3 assembly as in B. (D) Expression of soft and hard HGT candidates, and all other genes, in the nHd.2.3 assembly. A high-resolution version of this figure is available in SI Appendix.

Presence of fHGT candidate transcripts in poly(A)-selected RNA is strong evidence for eukaryotic-like expression and integration into a host genome. We mapped H. dujardini mRNA-Seq data (42) to UNC (Fig. 2C). Only nine of the UNC scaffolds that had low or no read coverage in our raw genome data had appreciable levels of mRNA-Seq reads mapped [between 0.19 and 31 transcripts per million (tpm)]. One of these (scaffold1161) contained two genes for which expression was >0.1 tpm, but all genes on this scaffold had best matches to Bacteria. The mRNA-Seq data thus gave no support for eukaryote-like gene expression from the low coverage, bacterial contigs in UNC.

Boothby et al. (13) assessed foreignness using an HGT index (47) and by analyzing phylogenies of candidate genes. However, these tests are only valid when there is independent evidence of incorporation into a host genome. Boothby et al. (13) assessed genomic integration of 107 candidate loci directly, using PCR amplification of predicted junction fragments. Most were adjudged confirmed, but no sequencing to confirm the expected amplicon sequence was reported. The 107 candidates were reported (13) to include 38 bacterial–bacterial and 8 archaeal–bacterial junctions (SI Appendix, Table S3). Our analyses identified 49 bacterial–bacterial junctions in their set, but these do not confirm HGT, as similar linkage would also be found in bacterial genomes. We found no expression of the 49 bacterial–bacterial junction loci, supporting assignment as contaminants rather than examples of fHGT.

Of the remaining 58 loci, only 51 were likely to be informative for HGT (SI Appendix, Table S3), as 7 were themselves or were paired with loci of unassigned taxonomic affinity. The informative loci included 24 prokaryotic to eukaryotic, 21 nonmetazoan eukaryotic–metazoan, and 6 viral–eukaryotic junction pairs. The UNC PacBio data confirmed only 25 of these junctions. All 58 loci had read coverage in Edinburgh raw data, and the same genomic environment was observed in nHd.2.3 for 51 loci (43 of which were HGT informative). We found evidence of expression from 49 of these loci.

The UNC H. dujardini genome is thus poorly assembled and highly contaminated. Scaffolds identified as likely bacterial contaminants in UNC included 9,872 protein predictions (Table 1). Evidence for extensive fHGT is absent, and most candidates were not confirmed by PacBio data, our read data, or gene expression. We present a more detailed examination of each of Boothby et al.'s claims for fHGT, including apparent congruence of codon use and presence of introns, in SI Appendix.

Low Levels of Functional Horizontal Gene Transfer in H. dujardini.

We screened nHd.2.3 for loci potentially arising through HGT. As mapping of transcriptome data to nHd.2.3 was equivalent to the precleaning nHd.1.0 assembly and better than the UNC assembly (SI Appendix, Table S1), the assembly has not been overcleaned. Forty-eight nHd.2.3 scaffolds (spanning 0.23 Mb and including 41 protein-coding genes) had minimal coverage in UNC data (Fig. 3C), suggesting that these were contaminants. The remaining 13,154 scaffolds spanned 134.7 Mb. Of the 23,021 protein coding genes predicted, only 13,500 had sequence similarity matches to other organisms, and of these, 10,161 had unequivocal signatures of being metazoan, with best matches to phyla including Arthropoda, Nematoda, Mollusca, Annelida, Chordata, and Cnidaria. A priori these might be considered candidates for metazoan–metazoan fHGT. However, as H. dujardini is the first tardigrade sequenced, this pattern may just reflect the lack of sequence from close relatives. The remainder had best matches in a wide range of nonmetazoan eukaryotes, frequently with metazoan matches with similar scores, and, for a few, bacterial matches. Some nonmetazoan eukaryote-like proteins (e.g., the two bodonid-like proteins discussed above) may have derived from remaining contamination.

We found 571 bacterial–metazoan HGT candidates in nHd.2.3, of which 355 were on 166 scaffolds that contained only other bacterial genes. Although some of these scaffolds also contained genes that had equivocal similarities, we regard them as likely remaining contaminants. Expression of these genes was in general very low (Fig. 3D), and we propose that these are “soft” candidates for fHGT. The remaining 216 HGT candidates were linked to genes with eukaryotic or metazoan classification on 162 scaffolds that had GC% and coverage similar to the tardigrade genome in both datasets (Table 2 and SI Appendix). Most of these (196, 0.9% of all genes) had expression >0.1 tpm (Fig. 3D) and are an upper bound of “hard” candidates for fHGT. However, phylogenetic analyses identified only 55 (0.2% of all genes) with bacterial affinities (having only bacterial and no metazoan homologs, or where analysis of alignments including the closest metazoan homologs confirmed bacterial affinities; SI Appendix). We identified 385 loci (1.7% of all genes) most similar to homologs from nonmetazoan eukaryotes (SI Appendix). Most (369) of these had expression >0.1 tpm, but phylogenetic analysis affirmed likely nonmetazoan origin of only 49 of these (0.2% of all genes; SI Appendix).

Within the high-coverage blob of assembly scaffolds supported by both Edinburgh and UNC raw data, blobtools analyses assigned 327 scaffolds as bacterial (black points in Fig. 3C). Fifty-two of these scaffolds were short (spanning 60.5 kb in total) and contained no predicted protein-coding genes, and 77 contained only predictions that were classified as eukaryote or unassigned. They were initially assigned as bacterial based on marginal nucleotide similarities to bacterial sequences. Many of the remaining 198 scaffolds were flagged in the gene-based analyses as containing fHGT candidates. Our assembly thus still contained contaminating sequences, mainly from bacteria but also including some from nonmetazoan eukaryotes. De novo joint assembly of the Edinburgh and UNC datasets in the future will permit robust elimination of such “difficult” contamination, as well as definition of the correct genome span, true gene content, and the contribution of HGT in H. dujardini.