Tracing RNA viruses associated with Nudibranchia gastropods

Introduction

With an estimated 10⁸ viral particles per millilitre of coastal water, oceans can be referred as hubs of viral diversity (Suttle, 2007; Mushegian, 2020). Here, viruses of unicellular organisms can govern the fate of host populations and influence global environmental cycles (Suttle, 2007; Brum et al., 2015; Sunagawa et al., 2020). As viruses are obligate intracellular parasites, their diversity reflects host genetic variety within the three domains of life (Krehenwinkel, Pomerantz & Prost, 2019). In this context, host-virus combinations frequently show co-adaptation processes, although viruses infecting multiple hosts and host-jumps have often been reported for invertebrate viruses (Shi et al., 2016). Completely leveraging the high-throughput power of modern sequencers, shotgun metagenomics is contributing to untangle the fascinating complexity of host and virus diversity and co-associations in the marine environment (Edwards & Rohwer, 2005; Carlos, Castro & Ottoboni, 2014; Kodzius & Gojobori, 2015). Starting from seawater filtrates or sediments subjected to RNA/DNA purification, amplification and subsequent high-throughput sequencing (HTS), viral metagenomics aims to reconstruct metagenome assembled genomes (MAGs) (von Meijenfeldt et al., 2019), often showing low similarity with database entries (Sullivan, 2015; Shi et al., 2016; Sunagawa et al., 2020). Although is itself challenging, the identification of viruses from metagenomic data is only an introduction to the understanding of host-to-virus associations and, finally, to the characterization of the functional roles of viral communities, also called viromes (Sullivan, 2015; Mara et al., 2020). Hence, both the distribution and the ecosystemic importance of viruses are probably underestimated and, to some extent, undiscovered (Koonin & Dolja, 2013; Garmaeva et al., 2019).

Differently from environmental metagenomics, the study of organism-associated biomes (holobiomes) can provide an overview of organism–organism associations (Hardoim et al., 2021) together with a snapshot of local viromes, as shown for the human virome (Garmaeva et al., 2019; Koonin, Dolja & Krupovic, 2021). Perhaps this approach is suffering less from limitations relating to proper host assignation, although additional challenges should be considered. Primarily, the low abundance of viruses within a healthy organism, together with the small size of viral compared to host genomes, makes their identification puzzling (Roux et al., 2015). Notably, some marine organisms can act as a magnifying glass for local viromes, due to their water filtration habits (Rosani et al., 2019). This has been exploited to reconstruct mimivirus genomes from oysters, as well as tracing human faecal footprints from mussel samples (Andrade et al., 2015; Olalemi et al., 2016). Although challenging, untangling host-associated viromes can also provide hints to trace the distribution of pathogenic viruses and to study host antiviral responses in non-model species (Rosani & Gerdol, 2017; Waldron, Stone & Obbard, 2018; Rosani et al., 2019).

Despite genomics, and with it metagenomics, having come on leaps and bounds over the past decades, several marine taxonomic groups have been overlooked resulting in a partial representation of the biological diversity of the sea (Greninger, 2018). One example pertains to the taxonomic group of gastropod nudibranchs, which includes more than 3,000 species seldom represented in sequence databases. These animals are characterized by colourful and varied body shapes, with sizes in the millimetric to centimetric range (Dean & Prinsep, 2017; Goodheart et al., 2018). Nudibranchs are known to derive defensive metabolites from the sponges they eat (Cheney et al., 2016), with several nudibranch-associated compounds that have been functionally studied highlighting their bioactivity against microbes and even viruses (Mudianta et al., 2016; Dean & Prinsep, 2017; Kristiana et al., 2020; Avila, 2020). Noticeably, viral-or bacterial-mediated diseases are not known for nudibranchs, whereas a fungus-like protist of the phylum Labyrinthulomycota has been associated with yellow and brown spot diseases observed in Tritonia diomedea (Collier et al., 2017). Infected nudibranchs mount an effective response able to encapsulate the pathogen with a lower impact on the fitness of the organism (McLean & Porter, 1982). In the context of global changes, nudibranchs have been proposes as sentinels of oceans. In particular warm-adapted nudibranch species appeared more susceptible to thermal changes, making these species highly vulnerable in the current temperature increasing scenario (Armstrong, Tanner & Stillman, 2019). RNA sequencing (RNA-seq) data are rare for these species, and most of them have been used to study the nervous system (Senatore, Edirisinghe & Katz, 2015), to trace nudibranch–gastropod phylogeny (Goodheart et al., 2017) or for comparative morphology studies (Goodheart et al., 2018). Noticing the complete absence of information regarding nudibranch-associated viruses, I exploited the transcriptomic datasets of 35 nudibranch species and, by performing a meta-transcriptomic analysis, I applied different pipelines to identify viruses associated with these species.

Materials and Methods

Results

Identification of the nudibranch-associated viruses

To identify bona fide viruses associated with these nudibranch RNA samples, DIAMOND blastx results were implemented with two different approaches: either used to identify the footprint of the viral RdRPs (palmscan) or to identify viral hallmark genes based on curated viral HMM databases. A total of 37 high-confidence RdRP footprints were identified by palmscan, whereas 128 putative viral contigs were identified by VirSorter2 (Table S1). Including the contigs identified by DIAMOND, a total of 552 contigs were classified as possible viruses, with only 19 contigs commonly identified by all the three methods (Fig. 1A, Data S1). The size distribution of the viral contigs increased from blastx to Virsorter2 results with the 19 common contigs placed in the upper part of the distribution (Figs. 1B and 1C). A taxonomic annotation of the 552 contigs based on the predicted CDS was performed with CAT, revealing that palmscan identified only contigs classified as ‘RNA virues’, whereas Virsorter2 identified also hits similar to metazoan genomes (Table S2). Like palmscan, Virsorter2 identified RNA viruses with the single exception represented by a Cossaviridae hit (ssDNA virus) identified due to an ORF encoding a parvovirus non-structural protein NS1 domain (SRR1950949_18608 in Data S1). Overall, several hits matching the phyla Kitrinoviricota, Negarnaviricota and Pisuviricota were identified, together with abundant unclassified viruses (Fig. 2 and Table S2). Regarding the 19 shared contigs plus the one identified by palmscan and Virsorter2, four contigs showed similarity to Beihai picorna-like virus 75 (Table 2). These contigs resembled three partial and one complete viral genome and their similarity to Beihai picorna-like virus 75 (KX883381) ranged from 46% to 53%. Two other contigs showed similarity to Wenzhou picorna-like virus 46, whereas four contigs were similar to different viruses belonging to the Pisuviricota phylum (Table 2).

Table 2:

Annotation of 20 viruses associated with nudibranch samples.

Contig ID	Size (bp)	NCBI nt database (DIAMOND)	Phylogenetic classification (CAT)				Reported host	Suggested host (palmID)
		Description	Phylum	Order	Family	Species
SRR3726693_23	6,219	Riboviria sp. isolate 3nj-RDRP-3				Beihai picorna-like virus 75	Echinoderm	Metagenome
SRR1950953_7108	9,998	Riboviria sp. isolate 3nj-RDRP-3				Beihai picorna-like virus 75	Echinoderm	Metagenome
SRR1950940_818	8,834	Beihai picorna-like virus 72	Pisuviricota			Beihai picorna-like virus 72	Anthozoa	Bacteria (Fig. 3A)
SRR1950953_2265	13,015					Beihai charybdis crab virus 1	Crustacea	Gastropod (Fig. 3B)
SRR1950953_444	10,163	Beihai picorna-like virus 75				Beihai picorna-like virus 75	Echinoderm	Metagenome
SRR3726693_131	3,162	Riboviria sp. isolate 3nj-RDRP-3				Beihai picorna-like virus 75	Echinoderm	Metagenome
SRR1950953_981	11,238					Wenzhou picorna-like virus 46	Pomacea canaliculata	Metagenome (viral - porifera)
SRR1950949_2406	9,654							Bacteria (Fig. 3C)
SRR3726693_286	10,990					Wuhan snail virus 2	Mollusca
SRR3726693_3931	9,200					Wenzhou picorna-like virus 46	Pomacea canaliculata	Metagenome (viral - porifera)
SRR3726693_6245	6,102	Beihai picorna-like virus 56	Pisuviricota			Beihai picorna-like virus 56		n.d.
SRR1950942_4074	8,133	Wenzhou picorna-like virus 19	Pisuviricota			Wenzhou picorna-like virus 19	Gastropoda	Plant (Fig. 3D)
SRR1950951_5889	6,244	Basavirus sp. isolate BaV/21164						Bacteria (Fig. 3E)
SRR3726705_1238	3,321							Bacteria (Fig. 3F)
SRR3726699_38607	4,326	Eptesicus fuscus rhabdovirus
SRR1950948_28126	3,196	Barns Ness breadcrumb sponge hepe-like 1	Kitrinoviricota	Hepelivirales	Hepeviridae	Barns Ness hepe-like virus 1	Porifera	Arthropoda (Fig. 3H)
SRR1950948_38692	1,519	Barns Ness breadcrumb sponge hepe-like 2	Kitrinoviricota	Hepelivirales	Hepeviridae	Barns Ness hepe-like virus 2	Porifera
SRR3726694_24261	3,010		Negarnaviricota	Mononegavirales	Nyamiviridae
SRR3726696_22111	4,656	Biomphalaria virus 4	Pisuviricota	Picornavirales		Biomphalaria virus 4	Biomphalaria	Bacteria (Fig. 3H)
SRR3726698_7138	1,701	Hubei earwig virus 1					Dermaptera

DOI: 10.7717/peerj.13410/table-2

Note:

For each contig, identified by the contig ID and size, the annotation obtained using blastn against the NCBI nt database (description) and the taxonomic classification obtained using CAT (phylum, order, family and species) were reported together with the suggested host of the best hit or by palmID. The viral contigs were ordered by total number of mapped reads. Underlined hits referred to genomes including low-confidence RdRPs (no palmID analysis was performed). Bolded taxonomic classifications referred to information not retrieved by CAT but added by palmID. SRR3726698_7138 represented the hit identified by palmscan and Virsorter2 only (annotations was retrieved using blastx).

Nudibranch-associated RNA virus distribution among the ‘planetary RNA virome’

A total of 15 high-confidence RdRP proteins were used to trace the distribution of the corresponding viruses within the collection of SRA datasets by using palmID implemented into the Serratus open science platform (Edgar et al., 2022). Association analysis based on the STAT kmer analysis performed on the RNA-seq datasets including similar viruses was used to retrieve possible host-virus associations (Fig. 3). The most abundant nudibranch-associated virus, Beihai picorna-like virus 75 (SRR3726693_23), was identified in metagenome-derived datasets poorly supporting possible hosts (Table 2). SRR1950953_226, which showed low similarity with a virus associated with crustacea (Beihai crab virus 1), was associated with arthropods by palmID. Differently, Wenzhou picorna-like virus 19 (SRR1950942_4074), initially associated with gastropods, is likely a plant virus. Viruses similar to Barns Ness hepe-like virus 1 and Barns Ness hepe-like virus 2, initially identified as possible porifera-associated viruses (Waldron, Stone & Obbard, 2018) but also found associated with bivalves (Rosani et al., 2019), revealed a possible association with arthropods. Beihai picorna-like virus 72 (SRR1950940_818) is possibly associated with Bacteria, as Biomphalaria virus 4 (SRR3726696_22111) and the unclassified SRR1950951_5889, SRR3726705_1238 and SRR1950949_2406 viruses. Notably, the two viruses similar to Wenzhou picorna-like virus 46 (SRR1950953_981 and SRR3726693_3931) were associated with a viral metagenome sample obtained from porifera (Rhopaloeides odorabile) as well as to several Pomacea canaliculata samples.

To further investigate possible host–virus associations in specific samples, a taxonomic classification of the assembled contigs of the eight samples including more than 500 viral reads was performed separately (Fig. S1). Except for nudibranchs, poor evidence of possible hosts of shared viruses was detected in independent samples (e.g. SRR1950953 and SRR3726693), nor for the highly abundant viruses found in sample SRR3726693, although bacteria are present in all these samples and might represent the genuine viral hosts (Fig. S1). An evident contamination of Dinophyceae (i.e. Symbiodinium spp.) is detectable in samples SRR1950951 and SRR3726699, but no viruses are shared between these samples. Low-level contaminations of Anthozoa are detectable in all samples, although the low abundance rejects the possibility that these species can act as hosts of the most abundant viruses. Differently, sample SRR1950953 showed Hydrozoa contamination, possibly representing the host of SRR1950953_226 and SRR1950953_981 viruses.

Discussion

Gastropod nudibranchs are intriguing models to investigate host–microbe associations and symbiosis from a functional point of view (Cheney et al., 2016; Mudianta et al., 2016; Kristiana et al., 2020). Notably, nudibranch-associated microbes are poorly known, and viruses have never being reported nor mentioned in the scientific literature. Besides the general interest of tracing viruses associated with these organisms, this will be preliminary to the understanding of possible roles mediated by viruses among nudibranch symbiosis as well as to untangle nudibranch antiviral defences. So far, no viral- or bacterial-mediated diseases are known in these species, possibly suggesting the presence of effective antiviral and antibacterial defences.

Here, the identification of viruses was approached by alternative and partially complimentary methods. Taxonomic classification of contigs (blastx) against the whole NCBI nr database identified both putative DNA and RNA viruses, although suffering from the poor representation of nudibranch species among databases. This aspect introduced abundant false positives in the identification of viruses, namely host contigs claimed as viral.

Differently, both palmscan and Virsorter2 genuinely identified RNA viruses. The former is designed to identify the footprint of the viral RdRPs (Edgar et al., 2022), an exclusive feature of RNA viruses. The latter identified viral hallmark genes and also provide a qualitative evaluation on the retrieved genomes (CheckV), reducing the identification of incomplete viruses, such as DNA viruses using RNA data. Overall, the sum of these methods identified 552 viral contigs, whereas the combination identified a core set of 19 viruses only. Sporadic DNA viruses identified in few samples and showing similarity with Herpesviruses, Mimiviruses or Poxviruses were discharged due to their genomic incompleteness.

Arguably, this result revealed the importance of using complimentary methods in order to detect high-confidence viruses associated with selected biological sample (e.g. nudibranchs), differentiating this analysis from a mere explorative survey focused on the detection of viruses. Likewise, although computationally very performant, the search of the RdRP only identified few viruses. Dedicated tools such as Virsroter2 limits false positives and increase the sensitivity by targeting multiple viral hallmark genes instead of the RdRP only, although requiring more computational effort.

By the search of the core nudibranch-virome among the ‘planetary RNA virome’ (Edgar et al., 2022), similar hits were retrieved as well as putative host-pathogen associations based on co-occurrence analysis.

Together, these results highlighted that nudibranchs possess different viromes, showing some similarities to viruses previously identified in gastropods, like Wenzhou picorna-like virus 46 and Wenzhou picorna-like virus 19 (Shi et al., 2016). Interestingly, the most represented virus showed similarity to Beihai picorna-like virus 75, initially found associated to echinoderms (Shi et al., 2016), although also present in several metagenomic datasets. Massive coverage in a single sample, together with the absence of evident signs of contamination of other possible host species, suggested that these viruses are hosted by nudibranchs (Learchis evelinae), and their role might be take into consideration in future research. Differently, all the other viruses are characterized by a lower coverage and co-occurrences analysis revealed that they possibly resembled bacterial viruses, which are still poorly studied (Wolf et al., 2020).

Indeed, nudibranch RNA samples are characterized by ‘physiological contaminations’, probably due to their feeding habits (corals) or mutualistic associations, here depicted by contigs of Symbiodinium spp. (Burghardt & Wägele, 2014; Mies et al., 2017). In this context, the host of the SRR1950940_818 contig similar to Beihai picorna-like virus 72 might be a coral, confirming its initial attribution to Anthozoa (Shi et al., 2016), although palmID associated also this virus with bacteria. Finally, the possible hosts of Barns Ness hepe-like virus 1 and 2 was considered, since these viruses were reported associated with Porifera and bivalves (Waldron, Stone & Obbard, 2018; Rosani et al., 2019), whereas they cooccur with Arthropoda (palmID) and nudibranch too.

Conclusions

A total of 20 RNA viruses were identified associated with nudibranch samples, confirming that RNA-seq datasets are useful data to untangle sample-associated viromes, particularly for under-studied taxa. While the analysis of RNA samples can reveal the active components of the biome as well as RNA viruses, DNA-based metagenomics will be required to trace the dormant DNA viruses. Starting from the evidence of RNA viruses associated with nudibranchs, future research can deepen the functional roles of viromes in these fascinating small and colourful marine species by untangling individual viromes using dedicated sample preparation approaches.

Supplemental Information