Analysis of the RNA virome of basal hexapods

Discovery of RNA viruses in public transcriptomic databases

Sequence database searches were finished in May 2019. The protein queries were RNA dependent RNA polymerase (RdRp) sequences representing every RNA virus family recognized by ICTV (Lefkowitz et al., 2018), as well as the majority of the RNA viruses that are unclassified. The protein queries were also sequences of structural proteins from diverse RNA virus families. The database analysed was the Transcriptome Shotgun Assembly (TSA) at the National Center for Biotechnology Information (http://www.ncbi.nlm.nih.gov). To detect all available representatives of the particular RNA viral family, database searches were performed iteratively. Comparisons were made using the TBLASTN program (Gertz et al., 2006), with the E-value cutoff set to 10⁻⁵ and default settings for other parameters. The most divergent representatives of the particular RNA viral family were used as queries. All newly obtained sequences were compared to reference protein sequences of all RNA viruses. Sequences yielding e-values larger than 1e⁻⁵ were retained and compared to entire NCBI NR database to exclude non-viral sequences. Sequences for which the top hit was a virus and sequences with no other BLASTP hits in NCBI NR Db were then treated as putatively viral in origin and subject to further analysis. To detect highly divergent viruses, we performed domain-based BLAST by comparing the newly obtained sequences against the conserved domain database with an expected value threshold of 1 × 10⁻². Sequences with positive hits to the RdRp domain were retained. DNA sequences were translated with the Translate program (web.expasy.org/translate/). The nucleotide sequences of all basal hexapod RNA viruses are available in the Data S1 file.

Analysis of endogenous virus elements

Endogenous copies of the RNA viruses were detected using the TBLASTN algorithm against basal hexapod genomes available in the Whole Genome Shotgun Database (WGS) at the NCBI, using viral protein sequences as queries. The queries involved protein sequences translated from both the virus genomes that were identified for the first time here as well as the reference virus genomes. Comparisons were made using the TBLASTN program (Gertz et al., 2006), with the E-value cutoff set to 10⁻⁵ and default settings for other parameters. For each potential endogenous virus, the query process was reversed to determine their corresponding phylogenetic group. The nucleotide and amino acid sequences of EVEs are available in the Data S1 file.

Prediction of protein domains

In order to recognize potential protein domains in the protein sequences analysed, we used NCBI CDD database (http://www.ncbi.nlm.nih.gov/Structure/cdd/wrpsb.cgi), by applying a cut-off E-value of 0.01. Some proteins were compared against SMART (smart.embl-heidelberg.de), InterPro (http://www.ebi.ac.uk/interpro/) and Pfam (pfam.xfam.org) protein domain databases at default parameters.

Phylogenetic analysis

To infer the phylogenetic relationships among RNA viruses, we used their RdRp protein sequences. Key representatives of the particular RNA viral family were included in the phylogenetic analysis. The protein sequences of the palm subdomain of RdRps were aligned using MAFFT (Katoh & Standley, 2013). Phylogenetic trees were reconstructed using the maximum likelihood (ML) method. For phylogenetic reconstruction, we used IQ-TREE with the in-built automated test to choose the best substitution model for each tree (Trifinopoulos et al., 2016). Branch support was computed for all trees using 100 replicates of parametric bootstrap, and 1,000 replicates of the approximate likelihood ratio test and ultrafast bootstrap. The iTOL online tool (http://itol.embl.de/) was used for phylogenetic tree annotation (Letunic & Bork, 2016).

Discovery of novel and highly divergent RNA viruses in basal hexapods

The collection of the NCBI Transcriptome Shotgun Assembly (TSA) and Whole Genome Sequence (WGS) databases for basal hexapods offers an attractive possibility to obtain the first insight into the diversity of their RNA viromes. We performed the analysis of RNA viruses in 16 transcriptomes and 6 genomes of the springtails (Collembola), silverfish (Zygentoma), diplurans (Diplura) and bristletails (Monocondylia: Archeaognatha) (Table 1). 12 transcriptomes and 2 genomes were positive for RNA viruses, while 4 transcriptomes and 4 genomes were negative. These data allowed us to identify ∼40 novel and diverse virus genomes or genome fragments that contained an RdRp domain. We observed extensive sequence divergence of the novel RdRp domains, most sharing 25–40% amino acid identity with previously described RNA viruses (Table 2).

The most complete set of the RNA viruses was obtained from the springtails, due to the largest number of the transcriptomes. All novel RNA viruses were compared with the known viruses in the NCBI databases. In such a way, we obtained information about the RNA virus family/clade they belong to and their similarity to the already known viruses. To infer their phylogenetic position in the particular RNA viral clade, we used ML phylogenetic analysis. Phylogenies of basal hexapod RNA viruses (Figs. 1–3) demonstrated that they belong to numerous RNA viral clades. The RNA viral genomes of basal hexapods have genome organizations that are very similar or differ only slightly from the winged insect representatives (Figs. S11–S16). Since the sample processing for the preparation of transcriptomic and genomic libraries of basal hexapods involved entire individuals, a substantial proportion of the viruses discovered here might be associated with undigested food, gut microflora or parasites that exist within the organisms investigated. However, homology searching and phylogenies show that basal hexapod associated RNA viruses are most closely related to the insect viruses. It should be noted that the RNA virus sequences identified in the analysed transcriptomes of basal hexapods are their “putative” viruses and specific experiments should be carried out to prove that these viruses are indeed replicating in these arthropod species.

Basal hexapods possess quite a diverse RNA virome

We found that basal hexapods possess representatives of 14 out of 24 RNA viral clades (Table 3). Such RNA virome diversity is higher than that of insect orders Blattodea (5/24), Dermaptera (6/24), Orthoptera (8/24), Lepidoptera (8/24) and Coleoptera (9/24). The only insect orders with similar or higher diversity of their RNA viromes are Odonata (12/24), Hemiptera (13/24) and Diptera (17/24) (Shi et al., 2016).

Basal hexapods possess three of the six known dsRNA viral clades: Reo, Partiti-Picobirna, and Toti-Chryso (Fig. 1). We found the first basal hexapod reovirus in Anurida maritima. Although this genome is partial, it is segmented. Six segments of reovirus were found and encode RdRp (VP1), VP2, VP3, VP4, VP5 and VP10 proteins (Data S1). This reovirus is quite divergent and shows less than 25% identity in the RdRp with the described reoviruses. It is most closely related to coltiviruses, significantly extending their host range (Fig. 1A). We found endogenized partitiviruses in the Machilis genome, they show 56% amino acid identity with the Culex mosquito partitivirus. Their RdRp has a surprisingly well conserved coding capacity. A partial sequence of the partitivirus RdRp was found in the transcriptome of Thermobia (Zygentoma); it has 37% amino acid identity with the Hubei partiti-like virus 10 (Fig. 1B). Toti-Chryso clade has a single representative in the transcriptomes of basal hexapods, in the Tetrodontophora springtail. We obtained only an RdRp fragment, and a few coat proteins. In the Toti-Chryso tree, the springtail representative groups together with the “diatom colony-associated dsRNA virus 10” (Fig. 1C). We found endogenized totiviruses in the Machilis genome.

Basal hexapods possess representatives of four clades of the negative-stranded RNA viruses: Mono-Chu, Bunya-Arena, Orthomyxo and Qinvirus (Fig. 2). In the Mono-Chu clade, we found only endogenized mononegaviruses in genomes of Machilis and Catajapyx that belong to chuviruses (Fig. 2A), nyamiviruses and rhabdoviruses (Data S1). An endogenized Bunyavirus nucleoprotein was found in the genome of Machilis (Monocondylia) and shows 20–30% amino acid identity with phleboviruses (Data S1). Especially interesting was the discovery of highly divergent representatives of orthomyxoviruses in Atelura and Catajapyx (Fig. 2B). Endogenous Catajapyx orthomyxovirus is represented only by the PB1 protein. These two novel orthomyxovirus PB1 proteins show just 30% identity with the known orthologs. In Atelura, we found a nearly complete orthomyxoviral genome encoding five of the six segments (PB1, PB2, PA, envelope and nucleoprotein) (Fig. S11). We found a full-length representative of the Qinvirus clade in springtails, in the Anurida. This sequence is quite divergent from the others reported recently (Shi et al., 2016), showing just 25% identity in the RdRp region. Anurida qinvirus extends the host range of this rare viral clade from a few protostomes to the basal hexapods (Fig. 2C). A complete genome of the Anurida qinvirus is 7,722 bp long, which represents a normal size for qinviruses. It is encoded in two segments, the larger one (5993 bp long) encodes RdRp, while the smaller one (1,729 bp long) encodes the putative structural protein that is homologous only to the Wuhan insect virus 15 (Fig. S12).

In basal hexapods, we found representatives of seven out of 12 positive-stranded RNA viral clades: Picorna-Calici, Hepe-Virga, Narna-Levi, Tombus-Noda, Flavi, Luteo-Sobemo and Permutotetra (Fig. 3). The largest diversity of the positive-stranded RNA viruses was found in the Picorna-Calici clade as the representatives of four picorna lineages were found –dicistrovirus (in Diplura only), iflavirus (in Diplura and Collembola), Kelp fly (in Zygentoma and Collembola) and Nora-like viruses (in Monocondylia only) (Fig. 3A). Few selected picornaviral genomes are shown in the Fig. S13. We discovered Hepe-Virga clade in Collembola and Diplura, where Negev-like viruses were prevailing. In the phylogenetic analysis of the Negevirus group, we included diverse representatives in basal hexapods, extending the host range of this RNA viral clade (Fig. 3B). The genome of the Sminthurus negev-like virus is quite large (∼10,6 kb). In addition to the RdRp-encoding ORF, it possesses additional ORFs with typical negevirus conserved protein domains (Fig. S14). We also found a complete Benji-like virus in Holacanthella (Collembola) transcriptome (Data S1). This virus is 45% identical in the core RdRp domain with the Hubei Beny-like virus 1, which was the first known metazoan benyi-like virus and was found only in Diptera (Shi et al., 2016). Its genome is 7,680 bp long and encodes a single ORF with 2,457 amino acids. The comparison of both metazoan benyi-like viruses showed that the springtail representative possesses a large region (between amino acids 850 and 1,969) that is absent in dipteran benyi-like virus. A single incomplete narnavirus genome (1,133 bp long) was found in springtails (Fig. 3C) and possesses a typical genome organization of narnaviruses (Fig. S14). Two representatives of the Tombus-Noda clade were found in Zygentoma, in Atelura and Tricholepidion transcriptomes (Fig. 3D). We found incomplete genomes of Tombus-like viruses in Zygentoma only (Fig. S16). In the Flavi clade, we found a number of fragments of jingmenvirus in a Sminthurus springtail, they encode both NS3 and NS5 proteins (Data S1). In a dipluran Megajapyx, we found few permutotetravirus fragments of the RdRp that show up to 59% amino acid identity with Hubei permutotetra-like virus 9 (Data S1). In the springtail Holacanthella, we found a fragment of sobemo-like virus RdRp that shows 47% amino acid identity with Hubei sobemo-like virus 17. In the transcriptome of the same species, we found sobemo-like capsid that shows 33% identity with Hubei sobemo-like virus 19. In the transcriptome of the Pogonognathellus springtail, we also found sobemo-like capsid (encodes viral-coat domain) that shows 34% identity with bat sobemovirus (Data S1).

Endogenous viral elements are not rare in basal hexapods

A considerable number of endogenous viral elements (EVEs) was discovered in insect genomes (Shi et al., 2016), which is an indication of past infection events (Holmes, 2011; Feschotte & Gilbert, 2012; Aiewsakun & Katzourakis, 2015). We searched for potential EVEs in all available basal hexapod genomes. We analysed RdRp proteins as well as numerous additional structural proteins, such as nucleo- and glyco-proteins, for representatives of all metazoan RNA viral families. EVEs in basal hexapods came from nine RNA viral clades: from Mono-Chu, Orthomyxo, Qin, Partiti-Picobirna, Reo, Tombus-Noda, Hepe-Virga, Bunya-Arena and Toti-Chryso (Figs. 1–3; Data S1). Since the RNA viromes of basal hexapods are diverse, it is interesting that they possess EVEs only in Diplura and Monocondylia genomes. Although arthropods possess EVEs for 13 RNA viral clades (Hepe-Virga, Luteo-Sobemo, Narna-Levi, Bunya-Arena, Mono-Chu, Orthomyxo, Nido, Partiti-Picobirna, Picorna-Calici, Reo, Tombus-Noda, Toti-Chryso and Qinvirus) (Shi et al., 2016), we found that the amount and the diversity of EVEs in basal hexapods is similar to them (Table 4). As expected, given their endogenous status, most of these sequences are only fragments of the parent virus genome (Figs. 1–3; Data S1). All EVEs in basal hexapods are integrated in random genomic loci in different species. In the vicinity of EVEs in basal hexapods no retroposon elements can be found. Our analysis demonstrated that EVEs are not rare in the genomes of some basal hexapods and have been generated by multiple independent integration events.

Composition and abundance of RNA viruses in viromes of basal hexapods

The comparison of basal hexapod RNA virome composition with that of insects (Shi et al., 2016) demonstrated some similarities and many differences. In both cases, Picorna-Calici clade is the largest. Hepe-Virga is the major additional clade in basal hexapods. The other RNA viral clades are much less abundant in the total RNA virome of basal hexapods, such as the Mono-Chu, Tombus-Noda, Narna-Levi, Partiti-Picobirna, Luteo-Sobemo, Orthomyxo, Reo, Toti-Chryso, Flavi, Permutotetra, Bunya-Arena and Qinvirus. Positive-stranded RNA viruses are prevailing in their RNA virome, while negative-stranded and double-stranded RNA viruses are less abundant in basal hexapods (Figs. 1–3). It was demonstrated that the abundance and composition of RNA viromes are obviously phylum-specific (Shi et al., 2016). While some RNA viral clades are much more abundant in winged insects, they are quite rare in basal hexapods. The reasons for such differences could be effects of biased sampling, depth and size of the RNASeq libraries, or real differences in the amount of some RNA viral clades. The majority of hexapods possess very similar patterns of RNA virome composition –few major RNA viral clades and numerous minor clades with limited distribution (Shi et al., 2016).

Comparison of winged insect (Pterygota) and basal hexapod RNA viromes

The basal position of apterygote hexapods in the hexapod tree (Misof et al., 2014) is important for understanding the origin and evolution of the insect-specific RNA viruses. We can compare novel basal hexapod RNA viruses with diverse relatives from winged insects. In such a way, we can trace the changes in the RNA viromes, originations of particular RNA viral families etc. (Table 5). It is obvious that basal hexapod and winged insect RNA viromes are similar, where insects collectively possess four viral clades more (18 of the 24 RNA viral clades in total) (Shi et al., 2016). It should be noted that all previously discovered insect RNA viruses were involved in the phylogenetic analysis of the invertebrate RNA virosphere (Shi et al., 2016).

dsRNA virome of the basal hexapods is represented by three RNA viral clades, while winged insects possess representatives of five viral clades. Some insect orders are without (Orthoptera), with a single (Lepidoptera, Dermaptera and Blattodea) or with just two dsRNA viral clades (Coleoptera). Insect orders with three or four RNA viral clades are Odonata, Hemiptera and Diptera. As evident from the Table 5, there are differences between the insect orders in the presence/absence of the particular dsRNA viral clade. A similar situation was observed in negative-stranded RNA viromes where basal hexapods possess representatives of four RNA viral clades. In insects, five negative-stranded RNA viral clades are present, but with unequal distribution patterns in diverse insect orders. Some of these RNA viral clades are diverse and rich (Mono-Chu and Bunya-Arena clades), while others are moderate (Orthomyxo) or very small (Ophio and Qinvirus). Some insect orders possess a single (Coleoptera), two (Lepidoptera and Dermaptera) or three RNA viral clades (Orthoptera, Blattodea, Odonata and Hemiptera). Diptera is the only insect order that possesses five out of six negative-stranded RNA viral clades. Until now, dipterans were the only insect order that possesses Quinvirus (Shi et al., 2016). Here, we show that qinviruses are indeed more widespread among insect orders (Fig. 2). The positive-stranded RNA virome in basal hexapods is represented by seven viral clades while winged insects possess eight clades out of twelve. In contrast to the basal hexapods, winged insects possess more diversity inside the RNA viral clades of the positive-stranded RNA virome. A number of insect orders possess a smaller number of positive-stranded RNA viral clades, such as Blattodea (1), Dermaptera (3), Lepidoptera (5), Orthoptera (5), Coleoptera (6) and Odonata (6). Hemiptera possess seven clades, while Diptera possess eight positive-stranded RNA viral clades. As in the double-stranded and negative-stranded RNA viromes, the distribution of positive-stranded RNA viral clades differs among insect orders. Only some of the RNA viral clades are widespread in hexapods, such as Picorna-Calici, Tombus-Noda, Mono-Chu, Bunya-Arena, Orthomyxo and Permutotetra. All other RNA viral clades have a much more limited distribution.

Horizontal virus transfer is very rare in basal hexapods

Many insects are known vectors for the dissemination of RNA viruses, such as mosquitoes and many plant pests (e.g., thrips, whiteflies, lepidopterans, coleopterans and scale insects) (Whitfield, Falk & Rotenberg, 2015; Rückert & Ebel, 2018). Since springtails (Collembola) are very abundant physical decomposers of plant and fungal material, there is a possibility of transfer of plant or fungal viruses into them. Springtails could potentially act as vectors of plant or fungal RNA viruses. However, the analysis of springtail transcriptomes and genomes showed that horizontal virus transfer (HVT) is extremely rare among them. We found only a single short fragment of a plant RNA virus in the springtail transcriptome, which could be present in ingested plant material infected with this virus. This was the alfalfa mosaic virus (AMV, Bromoviridae), found in the Holacanthella transcriptome (GFPE01073448, 340 bp long fragment, 99% amino acid identity to AMV). We were unable to find any sign of HVT in any other basal hexapod lineage.