Introduction

The first evidence that nematodes might host viruses were a number of reports of virus-like particles by electron microscopy (see Bird and Bird1 pp276–280 for a review of this early literature2). This work suggested that a number of free-living and plant parasitic nematode species could have viral infections3, but most early reports did not further characterise the viruses themselves, although Poinar and Hess4 identified a virus related to the iridoviruses in the cytoplasm of Romanomermis culcivorax. The recent discovery of three nematode nodaviruses5,6 from natural populations of Caenorhabditis elegans and C. briggsae, represented the first viruses known to naturally infect free-living nematodes and has rekindled interest in nematode virology.

Nodaviruses are positive sense single-stranded RNA viruses. This group was previously known to infect both vertebrate (fish) and invertebrate (arthropod) hosts, and the classical view of this family divides them into two clades (alpha- and beta-Nodaviridae) representing these divergent host groups. The nematode viruses thus significantly expanded the known host range of nodaviruses, and expanded the known diversity of the viruses themselves: the Caenorhabditis viruses form a group only distantly related to the arthropod nodaviruses. The isolation of viruses infecting the model nematode permits the investigation of nematode-virus interactions using a natural infection6,7. Research on these viruses has demonstrated that they infect intestinal cells in the worm7, have revealed a strikingly different capsid structure of Orsay (one of the C. elegans viruses) to the arthropod betanodaviruses8 and is starting to investigate the molecular biology9 and even develop approaches for genetic manipulation10 of this virus.

Following this discovery, a number of other potential nematode viruses have been described, significantly expanding the diversity of both nematode hosts and of the viruses infecting them. Four negative sense RNA viruses belonging to different families have been described in transcriptomic data from the plant parasitic nematode Heterodera glycines11. A survey of nematode EST data also revealed sequences related to the positive sense single-stranded RNA Picornavirales in two plant nematodes, Globodera pallida and Heterodera schachtii12. While the nature of these data, and while the presence of nematode-vectored plant viruses within the picornavirales13 make it impossible to be certain whether these are nematode viruses or contaminants from the plant material, it seems likely that these represent another novel group of nematode viruses.

While these findings have given us a glimpse of the possible diversity of nematode viruses, it seems likely that these rather sporadic and opportunistic findings may significantly underrepresent both the host range and diversity of viruses themselves: as far as we are aware, only two studies11,12 have systematically searched for virus-like sequence in a nematode, and these only examined transcriptomic data for RNA viruses. In other taxonomic groups, mining genomic data to find endogenous viral elements – DNA sequences derived from viruses following integration into the genome – has revealed a widespread and diverse menagerie of endogenised viral sequences14. These endogenous viral copies are, of course, particularly widespread for retroviruses where integration into the host genome is a obligate part of the lifecycle15 but integrated elements derive from a wide range of viruses, including viruses that have an exclusively RNA lifecycle. While the first report of endogenous elements derived from RNA viruses was only in 200416, there have subsequently been reports of at least 10 families in animal genomes14,17, including humans18 and other reports in fungi and bacteria (see Holmes et al.17 for review). Study of ‘fossil’ copies of viruses has given valuable insights into the evolution of both viruses19,20,21 and hosts22,23.

Here, we exploit the published sequence data for Caenorhabditis nodaviruses and published genomic data for nematodes to investigate the potential diversity of integrated nodavirus-derived elements in nematodes. We find an endogenous nodavirus-related sequence in the genome of Bursaphelenchus xylophilus24, the pine wood nematode: a plant parasitic nematode that causes severe damage to forestry and forest ecosystems25. The nodavirus-related sequence is embedded in a degerate LTR retrotransposon, suggesting a route via which this positive sense RNA virus could have become integrated into a nematode genome.

Results

Discovery of eBxnv-1

The search strategy we employed found only a single blastp hit with E-value < 10−2, in the genome of Bursaphelenchus xlyophilus24 to BUX.s01281.241 for all three RNA-dependent RNA polymerase sequences, with blastp E-values between 2 × 10−18 and 2 × 10−25 (bit scores between 91 and 113). The predicted gene model in B. xylophilus is a single exon gene encoding a hypothetical protein of 570 AA, with similarity to Prosite model PS5057 (RNA-dependent RNA polymerase of positive ssRNA viruses catalytic domain, at residues 407–529)26, and to the Superfamily HMM SSF56672 (DNA/RNA polymerases superfamily, at residues 407–529)27. A search using NCBI blast for sequence similar to this protein in the Genbank nr database found significant similarity only to RNA-dependent RNA polymerases of nodaviruses (Supplementary Table S3c). Tblastn identified only the same similarity regions as blastp, but blastn searches produced multiple small (around 60 bp) but significant hits for a region of Le Blanc virus to a number of Caenorhabditis species (4 hits in C. brenneri; 1 hit in C. briggsae; 61 hits in C. japonicum; 2 in C. remanei; 1 in C. sp 5 and 4 in C. sp 11). These sequences apparently form part of the internal transcribed spacer (ITS1) of the rRNA array in these species, and could represent capture of some host RNA by this virus, but the biological significance of these matches is unclear. Phylogenetic analysis confirms that the Bursaphelenchus predicted protein is related to nodavirus RDRP (Fig. 1). Intriguingly, there is strong support for this protein to be closely related to a subgroup of the arthropod-infecting alphanodavirus from Lepidoptera. The endogenous nodavirus we have identified – provisionally called eBxnv-1 – is thus only distantly related to the clade of Caenorhabditis nodaviruses.

Figure 1
figure 1

Phylogeny of eBxnv-1 with RNA-dependent RNA polymerase genes of other nodaviridae.

Values next to branches are bootstrap proportions for the partition implied by that branch, based on 1000 replicates. Bootstrap values for relationships within the beta-nodaviruses and between ANV (Drosophila melanogaster Amercan nodavirus), Flock house and Black beetle viruses are not shown.

Full size image

Structure of eBxnv-1

The gene model BUX.s01281.241 is significantly shorter than the polymerase proteins for other nodaviruses, at only 570 amino acid residues versus between 950 and 1050 for known virus proteins. In general, the amino acid sequence of these viruses is rather poorly conserved, and the eBxnv-1 is particularly divergent. Intriguingly, however, three conserved domains characteristic of nucleic acid polymerases28 and conserved in virus RNA-dependent RNA polymerases29 are clearly present in eBxnv-1 and align well against those in other nodavirus sequecnes (Supplementary Fig. S1; see PROSITE PS50507). In particular, the putatively catalytic Asp residues in domains IV and VI and Asn residue in domain V are present.

With any draft genome assembly, individual scaffolds could represent contamination rather than genuine sequence from the target organism. To confirm that the identified gene model is part of the B. xylophilus genome, we checked that the adjacent gene models show similarity to nematode proteins. While both predicted proteins 5′ to the eBxnv-1 locus show significant similarity to proteins from other nematodes, none of the three 3′ adjacent proteins have convincing BLAST hits to any proteins in the NCBI database (Supplementary Table S3a,b,d–f). Furthermore, a gap region between BUX.s01281.240 and BUX.s01281.241 makes the co-linearity of these genes in the genome assembly uncertain. We performed genomic PCR to fill this gap, with direct sequencing of the PCR products revealing the gap was of 1,828 bp. This genomic PCR confirms the correct assembly of this region, linking the identified endogenous sequence to these adjacent genes and to the rest of the scaffold (Fig. 2). Further investigation of the scaffold after filling this gap revealed the region has long terminal repeat (LTR) transposon-like structure (Fig. 2) with two LTRs, polypurine tract, primer binding site and pol coding sequences including peptidase and integrase between the two LTRs, but no sign of the gag proteins found in most full-length LTR transposons. This LTR appears to be a member of the Bel/Pao retroelements, and particularly to be related to members of the Tas clade of this group (Supplementary Table S3g,h), which includes elements described from both cnidarians and nematodes including Caenorhabditis elegans, Brugia malayi and Trichinella spiralis30. The predicted protein sequence show similarity to both the ribonuclease H and reverse transcriptase domains of the Pol open reading frame of these other elements, but are either only distantly related or degenerate. The LTR sequences themselves show complete conservation over 132 nucleotides and no homology outside this region – the Tas element has LTRs of around 250 nt, and the LTRs of all other elements in this clade are all longer. The RdRp-like protein is located within this LTR structure, along with an additional protein of unknown function (Fig. 2).

Figure 2
figure 2

(A)Structure of the putative endogenised nodavirus RNA1 locus in the Bursaphelenchus xylophilus genome, the surrounding LTR retrotransposon and adjacent gene models. (B) Genomic PCR confirming assembly correctness as this locus, and so confirming endogenous nature of this virus-derived gene. Lower case alphabets (a–d) on each lane represent amplified genome regions shown in (A). M; 2-log ladder DNA size maker.

Full size image

Expression of eBxnv-1-containing element

To investigate expression, RNA-seq reads were mapped to the gap-filled reference genome. The 100 bp-paired RNA-seq were mapped throughout the region between the two LTRs, but mapping coverage is as low as ~20x (FPKM ~60) (Supplementary Fig. S2A) and no significant expression difference between developmental stages (propagative-stages and dauer-stage) was observed. De novo assembly of the RNA-seq using Trinity generated a single contig of 7088 bp that again covers the entire LTR region. In addition, RT-PCR using mixed-stage nematodes detected as long as ~2 kb and ~5 kb fragments (Supplementary Fig. S2B). These results suggest this region is expressed at a low level as a single long transcript. Stage specific expression was further examined using RT-qPCR. We observed similar expression levels in adult females, males, L2 and L3/L4 although the expression levels are low compared to the actin gene (approx. 1/1000 of the level) used as an endogenous control in the qPCR experiment (Supplementary Fig. S2C). No virus particles were observed in a nematode homogenate under the electron microscope (Supplementary Fig. S3).

Other repeat elements and other copies of eBxnv-1 in the Bursaphelenchus genomes

A genome-wide, structure-based de novo search for LTR retrotransposons identified 1,301 potential candidate elements, including the element around eBxnv-1. 88 of those showed at least one HMM hit to a selection of endogenous retrovirus-related protein domain models from Pfam and GyDB. Only four show a set of domains which could be considered complete for transposition (RT, protease, integrase and RNaseH). Comparison of the eBxnv-1-containing element with the other candidates resulted only in short and low-specificity hits to 11 predicted candidates, suggesting they are not additional copies of the element in question. Genome-wide similarity searches for the eBxnv-1-containing element only delivered 9 insertion fragments (some on short contigs) and a number of isolated LTRs.

Next, we examined presence/absence of the eBxnv-1 element in other strains of B. xylophilus and in its closely related species by mapping of gDNA Illumina reads and PCR amplification using the specific primers (Supplementary Table S1). The element was detected in four strains of B. xylophilus (S10-P3, S10-P9, T4 and the reference strain Ka4C1) at the same genome location, but not in two other B. xylophilus strains (OKD1-F7 and C14–5) or B. mucronatus (Fig. 3), suggesting eBxnv-1-containing element were endogenised after the speciation of B. xylophilus and B. mucronatus. No SNPs and indels were detected in the area of eBxnv-1 (a total of 1713 bp) in the four B. xylophilus strains although their genomes are shown to be highly diverged31.

Figure 3
figure 3

Endogenisation of eBxnv-1 in B. xylophilus.

A maximum likelihood phylogenic tree for six B. xylophilus strains and the related species B. mucronatus based on a total of 1121 SNP sites. B. mucronatus was used as an outgroup. Plus signs indicate the presence of the eBxnv-1 element in that strain was identified by genomic sequence reads and PCR, minus indicates failure to detect the element. The arrow indicates the inferred timing of the integration of eBxnv-1 in B. xylophilus.

Full size image

Discussion

Nodaviruses are small, non-enveloped +sense ssRNA viruses, with RNA1 encoding the RNA-dependent RNA-polymerase and protein B2 (ecoded from the 3′ end of RNA1), and RNA2 encoding the capsid protein precursor32,33. Our results suggest that the RNA1 transcript of a nematode-infectious nodavirus infecting the germ line of B. xylophilus has been reverse transcribed by a cellular reverse transcriptase activity, captured by non-homologous recombination by an LTR retroelement and inserted into the host genome. This is the first evidence of a nematode nodavirus outside Caenorhabditis spp., expanding the known host-range of nodaviridae; the ancestral virus our element descended from represents a distinct lineage from the Caenorhabditis viruses or any other nodavirus. The Orsay virus apparently causes damage to intestinal cells in C. elegans6, so if an active nodavirus infective Bursaphelenchus was isolated, and shown to be pathogenic, this could represent a possible avenue for biological control of pine wilt disease34.

Our results are – to the best of our knowledge - the first confirmation of an endogenous viral element in a nematode genome, although the virus-related transcripts previously identified11,12 could be from endogenous loci. This builds on the growing evidence of extensive endogenised virus elements found by in-silico screening of mammal, bird and insect genomes against a library of non-retroviral peptides using methodology similar to that employed here14. As commonly found in EVEs in other systems14,21, the endogenous copy we identify is not closely related to any of the exogenous viruses within this family, and the most closely related viruses use insect, not nematode, hosts. This finding thus expands our knowledge of the diversity of nodaviridae and of the diversity of host-virus relationships in this group. In common with other recently discovered EVEs related to ssRNA+ viruses, eBxnv-1 appears to be present in very low copy number – in this case, a single copy – perhaps due to low viral mRNA copy number17.

The entire LTR element – including the nodavirus-derived ORF – is expressed at a low level across all B. xylophilus life history stages investigated. Expression of the endogenised viral element, and the conservation of several key domains required for polymerase function in an otherwise divergent amino acid sequence, raises the possibility that this gene is functional within nematode cells. A different virus RdRp gene has recently been propose to be functioning in some bat genomes35 – while it is hard to be sure about possible functions, these authors speculate that this mammalian protein could play a role in either anti-viral defense or in the RNAi pathway. While rare, a number of examples are known where viral proteins have been ‘domesticated’ by the cellular host, particularly being recruited to anti-viral defence. The most well-explored examples are in mammals (see ref. 23 for review) but this process is likely to be more widespread (e.g. ref. 36). The fact that the only appearance of this LT element in the present-day B. xylophilus genome is the copy enclosing eBxnv-1 makes is tempting to speculate that there could be some functional explanation for the persistence of this locus, but we have no evidence for any such function.

The presence of a +sense RNA virus integrated into the Bursaphelenchus genome implies that these nodaviruses can infect germline cells and that genetic material derived from these viruses can find its way into the nucleus. As the nodavirus has no reverse transcriptase activity, the ancestral viral RNA of our eBxnv-1 element must have interacted with an exogenous protein with reverse transcriptase activity to generate a DNA molecule14,17,18. In this case, we presume that this RT was from the LTR element we now find surrounding the eBxnv-1 protein itself and speculate that the recombination between the retroelement and nodavirus templates occurred during the reverse transcription process itself by RT copy-choice recombination, and that this recombinant molecule could then integrate into the host DNA. This process of retroelement capture of a non-retrovirus material has been experimentally demonstrated in a mammalian system37, but we are aware of only one previous natural example38, and in that case only the terminal repeats of the implicated element remained.

The LTR retrotransposon we identify in the B. xylophilus genome is clearly derived from an element related to Bel/Pao elements previously identified in other nematodes. The long terminal repeat sequences themselves appear to be truncated but a well-conserved region is still present. However, because the protein-coding open reading frames appear to be highly degenerate, the LTR element is certainly non-functional. The element appears to be restricted to one clade of B. xylophilus strains and absent in the related species B. mucronatus, suggesting that the endogenisation event occurred relatively recently in B. xylophilus. It is noteworthy that the endogenised region in the element (i.e. eBxnv-1 gene) retains key catalytic sites and is highly conserved in the four B. xylophilus strains. Taken together with the fact that the region is expressed in multiple developmental stages of the nematode, albeit weakly, suggests that the RdRP gene could be functional although the evolutionary significance for the host is unclear. Similar expression and possible exaptation of endogenised RdRp has been observed in Eptesicus bats18.

Methods

BLAST searches (blastp, tblastn and blastn) with the nucleotide and predicted protein sequences for all 3 nodavirus genomes available against the genomes and predicted proteomes for all 15 nematode genomes available in WormBase release WS23339 (7 Caenorhabditis spp., Ascaris suum, Brugia malayi, Bursaphelenchus xylophilus, Haemonchus contortus, Meloidogyne hapla, Pristionchus pacificus, Strongyloides ratti and Trichinella spiralis, with genomic sequence only (no predicted proteins) available for Heterorhabditis bacteriophora and Meloidogyne incognita). Protein sequences for the RNA-dependent RNA polymerase protein of other known nodaviruses were downloaded from Genbank (see Supplementary Table S1), and aligned with mafft v6.85740 with the ‘—auto’ option. The alignment was then cleaned with trimAl v1.4 using the ‘-automated1’ option41. Phylogenetic analysis was performed using RAxML v.7.2.842, using a model of amino-acid substation (LG+F) chosen as that minimizing AICc and a discretized gamma distribution to model rate variation across sites plus an estimated proportion of invariant sites. Inference used the rapid hill-climbing algorithm of RAxML.

An assembly gap between BUX.gene.s01281.240 and BUX.gene.s01281.241 was filled by PCR direct sequencing. PCR amplification was performed using Tks Gflex DNA Polymerase (Takara, Japan) as recommended by the manufacturers; 30 cycles of PCR (in 30 μl volume) were performed with 10 sec at 98 °C and 7 min at 68 °C with 2 min of 94 °C initial denaturation. PCR products were sequenced on ABI Prism 3130 DNA sequencer using BigDye Terminator v3.1 (Applied Biosystems). The primers used for PCR and sequencing are listed in Supplemental Table S2. RNA-seq reads were obtained from poly-A selected RNA extracted from mixed-stage and dauer-stage B. xylophilus43 using Illumina HiSeq2000 platform and standard protocols (http://www.illumina.com/). Reads were mapped against the reference genome using TopHat software, v 2.0.6 (–mate-std-dev 50 -a 6 -i 10 -I 20000 –microexon-search –min-segment-intron 10–max-segment-intron 50000)44. De novo RNA-seq assembly was performed using the Trinity pipeline with default options45. For RT-qPCR, total RNA was extracted from each stage of nematodes using TRI reagent (Sigma). cDNA was synthesized using iScript according to manufactures instruction (BioRad). RT-qPCR was performed using 1 μl cDNA and Power SYBR Green PCR Master Mix (Life technology) on the StepOnePlus System (Applied Biosystems).

To search for virus-like particles by electron microscopy, approximately 100,000 B. xylophilus individuals (strain Ka4, the original reference genome strain which was isolated in Ibaraki, Japan in 1994 and has been maintained by subculturing in FFPRI) were homogenized in 0.1 M sodium phosphate buffer at pH7.5. Cell debris was pelleted by low-speed centrifugation (8,000 g). The supernatant was adjusted to 6% PEG6000 and 0.125 M NaCl. After centrifugation at 10,000 g, the pellet was suspended in 0.03 M sodium phosphate buffer, pH7.5 and subjected to ultracentrifugation at 100,000 g. Extracts were negatively stained with 2% phosphotungstic acid (pH6.5) and examined in a Hitachi JEM2000EX transmission electron microscope at an 80 kV acceleration voltage.

LTRharvest46 was used to initially detect the presence of LTRs flanking the virus insertion (LTR lengths 10–1000 bp, total length 1–15 kbp, 90% minimum LTR similarity, TSD length 5–20 bp, seed 100 bp, looking for motif starting ‘tg’ and ending in ‘ca’). To assess the abundance of copies of the detected LTR retrotransposon in the B. xylophilus genome, LTRharvest was run with more sensitive parameters on the whole genome sequence (LTR lengths 100–1000 bp, total length 1–10 kbp, 85% minimum LTR similarity, TSD length 4–20 bp, seed 30 bp) to identify potential candidates in the genome. LTRdigest47 was used to subsequently annotate the candidates with protein domains from Pfam48 and GyDB30 as well as putative primer binding sites and polypurine tracts. The additional information was processed using custom scripts to remove low confidence candidates by only keeping (1) candidates with at least one protein domain, and (2) candidates with full LTR retrotransposon domain sets (reverse transcriptase, protease, integrase, RNAse H). The sequence of the nodavirus containing LTR element was then compared to the remaining genome-wide LTRharvest results as well as against the whole genome sequence using high sensitivity (E-value: 0.01) BLASTN and TBLASTX searches. Classification of the LTR element was by BLAST similarity search of the LTRs and LTRharvest protein matches against GyDB release 2.0, using the tool available at http://gydb.org/index.php/Blast.

To examine presence/absence of the eBxnv-1-containing element, genomic DNA reads from six B. xylophilus strains (Ka4C1, T4, S10-P3, S10-P9, C14-5 and OKD-F7)31 and B. mucronatus (Un1 strain) was mapped to the B. xylophilus v1.2 assembly with the gap-filled scaffold01281 using SMALT v0.7.4 (https://www.sanger.ac.uk/resources/software/smalt/) with options –x and –y 0.8. Variants in the scaffolds were called using the Picard tool (ver. 1.95) and GATK (version ver. 3.3.0) as described previously31. A maximum-likelihood tree was generated from SNP sites derived from 10 biggest scaffolds in the assembly using SNPhylo program49 with options (-m 0.1 –M 0.1 –l 0.9 –b 100). Presence of the element was confirmed by PCR amplification using genomic DNA and primers listed in Supplementary Table S2.

Additional Information

Accession Codes: Nucleotide sequence of the eBxnv-1 region has been deposited to DDBJ/GenBank/EMBL under accession number LC158686. The authors have no competing interests to declare. B. xylophilus population sequence data, B. xylophilus RNA-seq data and DNA read data of B. mucronatus was deposited at DDBJ Sequence Read Archive (DRA) under BioProject numbers PRJDB3459, PRJDB3458 and PRJDB3460, respectively.

How to cite this article: Cotton, J. A. et al. An expressed, endogenous Nodavirus-like element captured by a retrotransposon in the genome of the plant parasitic nematode Bursaphelenchus xylophilus. Sci. Rep. 6, 39749; doi: 10.1038/srep39749 (2016).

Publisher's note: Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.