Diverse RNA Viruses Discovered in Three Parasitoid Wasps of the Rice Weevil Sitophilus oryzae

The enormous diversity of RNA viruses in insects is continuously validated. Parasitoid wasps, as biocontrol insects which are widely used against insect pests in agroecosystems, may also carry many “good” RNA viruses.

stranded RNA viruses (1)(2)(3). Species abundance and diversity of RNA viruses also exist in insects, the largest group in invertebrates (4)(5)(6). Nowadays, the application of metagenomics accelerates the studies of viral identification, persistence, spread, and interaction with both their insect hosts and other microbes. For example, using metagenomic analysis, seven novel positive-strand RNA viruses, belonging to the families Dicistroviridae, Parvoviridae, and Circoviridae, are identified in an important stingless bee Melipona quadrifasciata which have the largest copy numbers in unhealthy bees and may be pathogenic (7). Eight novel RNA viruses, including four positive-strand RNA viruses, two negative-strand RNA viruses, and two double-stranded RNA, are discovered in the high-throughput sequencing data of an important invasive agricultural insect pest, the oriental fruit fly Bactrocera dorsalis (8). Through a metatranscriptomics analysis, many new RNA virus species are discovered in different mosquito populations, including many that affect human health. This indicates that the diversity of RNA viruses is affected by the host mosquito taxon and genetic background (9)(10)(11). With understanding of mosquito viromes, the prevention and control of pathogenic arboviruses will be more easily manipulated.
Parasitoid wasps (Hymenoptera: Apocrita) are widely used as biocontrol agents against insect (Arthropoda: Hexapoda: Insecta) pests in agroecosystems. They are frequently associated with viruses or virion-like particles (VLPs). In the past decade, RNA viruses have been reported in some parasitoid wasps, as well as in virus-host interactions. Several viruses have been found to affect the host wasps or the wasps' host. For example, Pteromalus puparum negative-strand RNA virus 1, an artovirus found in Pteromalus puparum wasp which preys on butterflies of several species, increases the longevity of the wasp while it reduces female offspring numbers (12).

RESULTS
Virus-like contigs obtained from the transcriptomes of three parasitoid wasps. Through transcriptome analysis, 29 virus-like fragments were obtained in three weevil wasp species, including A. calandrae (8), L. distinguendus (3), and T. elegans (18) ( Table S3). The majority of small contigs were in T. elegans. Based on BLASTX results, these fragments were classified into one of two main groups: positive-strand RNA virus (11) or negative-strand RNA virus (18) (Table S3). All three wasps had these two types of virus-like sequences (Table 1 and Table S3).
WWPSRV-1, WWPSRV-2, and AcPSRV-1 were related to positive-strand RNA virus. (i) WWPSRV-1. Three large contigs obtained from the wasps' transcriptomes (AC_Contig_1, LD_Contig_1, and TE_Contig_1) (Text S1) were highly related to Nasonia vitripennis virus and could be assembled into one big virus-like contig (12,325 bp), named WWPSRV-1 (Table S3 and Fig. S1A). By sequencing the assembled viral genome, including the 59 and 39 genome termini, the complete genome of WWPSRV-1 was confirmed to be 12,332 nucleotides (nt) in length, excluding poly-A ( Fig. 1A and 2A). G1C pairs comprised 39.94% of the nucleotides. The WWPSRV-1 genome contained four ORFs in 59to 39 orientation, located at nt positions 1,003 to 8,739, 8,711 to 9,385, 9,367 to 10,947, and 10,940 to 12,007, respectively ( Table 2). The four deduced ORFs of WWPSRV-1 were not in the same frame. ORF1 and ORF3 were in the same frame, while ORF2 and ORF4 were in the other frame. ORF1 to 4 were overlapping by a fragment (range from 2 to 23 bp, excluding the start and stop codons), respectively ( Table 2 and Fig. 2A). The leader and trailer regions of the WWPSRV-1 genome were 1,002 and 325 nt in length, respectively (Table 2). A palindromic structure existed in the 59 leader regions ( Fig. 2A). The trailer regions ended with poly-A. Transcriptome sequencing (RNA-seq) reads of the three wasps mapped to the WWPSRV-1 genome exhibited similar fluctuating distribution on the viral genomic RNA, with a lower coverage close to 59 ends of the genome ( Fig. 2A).
(ii) WWPSRV-2. Four contigs (AC_Contig_3, LD_Contig_2, and TE_Contig_2 and 3) (Text S1) were assembled into a big virus-like contig (4,899 bp), which had a high similarity with Diabrotica virgifera virgifera virus 1 (DvvV-1) polyprotein (coverage range from amino acid [aa] 1,440 to aa 2,995) (Table S3 and Fig. S1B). Moreover, AC_Contig_2 (3,308 bp) (Text S1) was also aligned to DvvV-1 polyprotein (coverage range from aa 210 to aa 1,002) (Table S3). By sequencing these contigs, the five contigs were verified to belong to the same virus, named WWPSRV-2 (Fig. 1B). The complete genome of WWPSRV-2 was 9,409 nt in length, excluding poly-A. G1C pairs comprised 40.95% of the nucleotides. Only one large ORF was predicted in WWPSRV-2 genome, located at nt 34 to 9,228 in 59 to 39 direction. The 59 and 39 untranslated regions (UTR) of WWPSRV-2 genome were 33 and 181 nt in length, respectively ( Fig. 2A and Table 2). The fluctuating distributions of the three read maps were similar, with the higher coverage area near 39 ends of the viral genome. However, a small number RNA-seq reads of L. distinguendus could map to the WWPSRV-2 genome (Fig. 1A). Of note, PCR detection results revealed that this virus could exist in the adults of S. oryzae ( Fig. 1B and Fig. S2).
(iii) AcPSRV-1. AC_Contig_4 (Text S1) was 7,558 bp in length, including poly-A, and contained regions similar with viral RNA-dependent RNA polymerase (RdRp) of Dicistroviridae sp. (Table S3). It contained two large ORFs in the 59 to 39 direction, located at nt positions 183 to 4,343 and 4,660 to 7,302, respectively. The contig was provisionally named AcPSRV-1, and G1C pairs comprised 34.10% of the nucleotides. The 59UTR, intergenic region (IGR), and 39UTR of AcPSRV-1 genome were 182, 242, and 247 nt in length, respectively ( Fig. 2A and Table 2). However, this contig failed to be confirmed by PCR. Only the RNA-seq reads of A. calandrae could map to AcPSRV-1 genome. The higher coverage area also appeared close to 39 ends of the viral genome ( Fig. 2A).
AcNSRV-1, AcNSRV-2, and LdNSRV-1 were related to negative-strand RNA virus. (i) AcNSRV-1. Four contigs related to negative-strand RNA virus were discovered in A. calandrae. The two long contigs (AC_Contig_5, 12,699 bp; AC_Contig_6, 12,118 bp) (Text S1), which both had high identities with Linepithema humile rhabdolike virus 1, were provisionally designated AcNSRV-1 and AcNSRV-2, respectively (Table S3). By sequencing, the confirmed complete genome of AcNSRV-1 was  (Table 3 and Fig. 2B). The fluctuating distribution of RNA-seq reads of A. calandrae exhibited that reads in the five ORF regions were much more numerous than those in the nearby untranslated region. The position of the peak corresponded exactly to the position of each ORF (Fig. 2B). The leader and trailer regions of the AcNSRV-1 genome were 208 and 437 nt in length, respectively. Their terminal nucleotides were complementary (Table 3 and Fig. 2B) and could form a putative panhandle structure which may be involved in genome replication (17).
(ii) AcNSRV-2. The complete genome of AcNSRV-2 failed to be confirmed by PCR. The known length was 12,118 nucleotides and G1C pairs comprised 37.85% of the nucleotides. The genome contained six ORFs in the 39 to 59 direction, located at nt positions 29 to 1,330, 1,518 to 2,624, 2,803 to 3,183, 3,453 to 5,075, 5,214 to 5,429, and 5,572 to 11,886, respectively (Table 3 and Fig. 2B). The six ORFs were also confirmed by the RNA-seq reads fluctuating distribution (Fig. 2B). The leader and trailer regions of the AcNSRV-2 genome were 28 and 232 nt in length, respectively. Their terminal nucleotides were not complementary.
(iii) LdNSRV-1. LD_Contig_3 (10,832 bp) in L. distinguendus, designated LdNSRV-1, showed a certain sequence similarity with a rhabdovirus, Wuhan ant virus. The complete genome of LdNSRV-1 verified by sequencing was 11,633 nucleotides in length with G1C  the  Table S1. pairs comprising 45.64% of the nucleotides ( Fig. 1D and 2B). It contained five large ORFs in the 39 to 59 direction, located at nt positions 44 to 1,408, 1,499 to 2,587, 2,611 to 3,342, 3,542 to 5,080, and 5,169 to 11,519, respectively (Table 3 and Fig. 2B). The fluctuating distribution of RNA-seq reads of L. distinguendus also indicated the existence of the five ORFs (Fig. 2B). The leader and trailer regions of the LdNSRV-1 genome were 43 and 114 nt in length, respectively. Their terminal nucleotides were not complementary.
In T. elegans, the contigs related to negative-strand RNA virus were too small to get the complete virus genome (Table S3).
WWPSRV-1 represented a novel virus related to Nora virus. General properties of all ORFs in WWPSRV-1 genome were listed in Table 2. The predicted translation products of WWPSRV-1 ORF1 to 4 were found to have amino acid sequences similar to those of the Nasonia vitripennis virus RdRp and ORF3, 4, and 5, respectively. ORF1 protein was predicted to be a large transmembrane polyprotein (2,578 aa and isoelectric point [pI] of 8.05) with three noncytoplasmic domains, two cytoplasmic domains, and four transmembrane regions (Fig. S3A). A specified putative helicase domain (aa 824 to 1,005), a trypsin-like serine protease domain (aa 1,773 to 1,995), and a RdRp domain (aa 2,157 to 2,530) were successively existent in the noncytoplasmic domain of ORF1 protein ( Fig. 2A and Fig. S4).
To determine the relationship of WWPSRV-1 with other picorna-like viruses, we conducted a maximum likelihood phylogenetic analysis based on the amino acid core   sequences of RdRp. Notably, WWPSRV-1 clustered as a distinct lineage in Nora virusrelated clade, indicating the need for a novel virus genus or family (Fig. 3). WWPSRV-2 was a novel iflavirus. Only a big ORF was predicted in WWPSRV-2 genome, which encoded a large polyprotein. The putative domains in the polyprotein were arranged successively in the order CP1, CP2, CP3, helicase, protease, and RdRp, corresponding to that in iflaviruses (Fig. 2A).
The polyprotein (3,064 aa and pI = 6.19) was a large transmembrane polyprotein that consisted of four noncytoplasmic domains, three cytoplasmic domains, and six transmembrane regions. All three predicted picornavirus-like capsid protein domains appeared in the first noncytoplasmic domain. A typical helicase domain (aa 1,440 to 1,611) was present in the third noncytoplasmic domain. Moreover, a relatively small motif (aa 1,799 to 1,842) predicted as a zinc finger domain was also in that noncytoplasmic domain. The putative protease (aa 2,212 to 2,417) and RdRp (aa 2,490 to 2,931) domains were in the last noncytoplasmic domain ( Table 2 and Fig. S3B and S4).
To determine the phylogenetic relationship of WWPSRV-2 in other iflaviruses, phylogenetic analysis was performed using the highly conserved RdRp region. WWPSRV-2 could form a clade with the known members of iflaviruses, such as sacbrood virus (SBV) and Lygus lineolaris virus 1 (LyLV-1) (Fig. 3).
AcPSRV-1 was a novel dicistrovirus. Two big ORFs were present in the AcPSRV-1 genome, which encoded nonstructural and structural proteins, corresponding to that in dicistroviruses. The nonstructural protein encoded by ORF1 was a large transmembrane polyprotein (1,386 aa and pI = 8.04) that consisted of two noncytoplasmic domains, one cytoplasmic domain, and two transmembrane regions. The helicase (aa 67 to 248), protease (aa 595 to 831), and RdRp (aa 963 to 1,350) domains were in the different two noncytoplasmic domains ( Fig. S3C and S4). The structural protein (880 aa and pI = 6.59) encoded by ORF2 was a polyprotein in which three picornavirus-like capsid protein domains exist ( Fig. 2A and Table 2).
Phylogenetic analysis of the highly conserved RdRp region was performed to confirm the relationship of AcPSRV-1 in other dicistroviruses. AcPSRV-1 could be clustered with Himetobi P virus (HiPV), black queen cell virus (BQCV), and Triatoma virus (TrV), which belong to genus Triatovirus (Fig. 3).
Gene junction analysis of the three negative-strand RNA viruses. By comparing the 59 and 39 untranslated regions and intergenic regions of the AcNSRV-1, AcNSRV-2, and LdNSRV-1, we pinpointed for each ORF the putative termination signals (TS), intergenic spacer (IS), and transcription initiation (TI). A highly conserved motif was identified in the 59 and 39 untranslated regions and intergenic regions of the AcNSRV-1, AcNSRV-2, or LdNSRV-1, respectively.
In AcNSRV-2 genome, the conserved transcription initiation motif of 39-UUCA(G/C)-59 was identified upstream of each putative ORF, and the transcription termination motif of 39-UAAAUCUUUUU(U)-59 was detected downstream of every ORF. The conserved intergenic spacer 39-GAG-59 was found in the intergenic region among N/P, P/ M, and M/L (Fig. 4B).
AcNSRV-1 and AcNSRV-2 were novel viruses belonging to family Lispiviridae. The general properties of the tentative proteins encoded by the AcNSRV-1 and AcNSRV-2 genome were shown in Table 3. Both the ORF1 proteins encoded by the two viruses had close to 30% identities with the putative capsid of Linepithema humile rhabdo-like virus 1. Through the VIRALpro analysis, the ORF1 proteins were suggested to have the capsid sequences (distance was 0.5556 in AcNSRV-1 and distance was 1.2285 in AcNSRV-2). Accordingly, the ORF1 proteins of the two viruses were supposed to be the putative nucleocapsid proteins (N), which had similarity in length (1,290 aa and 1,302 aa, respectively) and isoelectric point (5.75 and 6.43, respectively).
Only ORF4 proteins encoded by the two viruses were predicted to be transmembrane polyproteins and to possess a cleavable signal peptide. Compared to the three transmembrane regions present in the AcNSRV-1 ORF4 protein, only one was in the AcNSRV-2 ORF4 protein ( Fig. S5A and B). The AcNSRV-1 ORF4 protein had 25.5% identity with the putative glycoprotein of Hubei rhabdo-like virus 3. In addition, it was classified to the fusion glycoprotein F0, Paramyxoviridae (E value = 1.00E208), within the InterProScan. The AcNSRV-2 ORF4 protein showed 25.12% identity with the putative glycoprotein of Blattodean arli-related virus OKIAV101. Moreover, only the ORF4 proteins of the two viruses were recommended to have the potential N-linked sites. All of these results indicated that ORF4 encoded the viral glycoprotein (G).
Phylogenetic analysis of the highly conserved region of L protein was performed to confirm the relationship of AcNSRV-1 and AcNSRV-2 in other lispiviruses. AcNSRV-1 and AcNSRV-2 could be clustered with Berant virus and Linepithema humile rhabdolike virus 1, which belong to the family Lispiviridae (Fig. 5).
LdNSRV-1 was a novel rhabdovirus. The general properties of the tentative proteins encoded by the LdNSRV-1 genome were shown in Table 3. The LdNSRV-1 ORF1 protein had a close-to-39.81% identity with nucleocapsid protein of hymenopteran rhabdo-related virus OKIAV109 (HyRRV-OKIAV109). It was suggested to have the capsid sequence (distance = 0.3497) with VIRALpro. Accordingly, the ORF1 protein was supposed to be the N protein. ORF4 protein was predicted to be a transmembrane polyprotein which only contained one transmembrane region and to possess a cleavable signal peptide (Fig. S5C). It has 28.07% identity with HyRRV-OKIAV109 glycoprotein. Moreover, ORF4 protein was recommended to have two potential N-linked sites. All of these indicated that ORF4 encoded the viral G protein. The LdNSRV-1 ORF5 was supposed to encode the viral L protein which also contained the putative RdRp, PRNTase, and MT domains (Fig. 2B and Fig. S6). It exhibited more than 40% similarity with HyRRV-OKIAV109 RdRp-complex. The LdNSRV-1 ORF2 protein (pI = 6.98), containing 40 putative phosphorylation sites, was most likely a P protein. ORF3 protein (pI = 8.94), containing 21 putative phosphorylation sites, may be an M protein.
Phylogenetic analysis of the highly conserved region of L protein was performed to confirm the relationship of LdNSRV-1 in other rhabdoviruses. LdNSRV-1 could be clustered with Lasius neglectus virus 2 and hymenopteran rhabdo-related virus OKIAV38, which belong to the genus Alphahymrhavirus (Fig. 5).

DISCUSSION
More and more RNA viruses have been discovered with the increasing transcriptome data of insects. Also, whole-transcriptome analysis has become a very important method in virus identification and surveillance in insects relevant to economics and public health. For example, 42 putative novel picornaviruses are discovered in varroafree Australian Apis mellifera (18), three novel RNA viruses are identified in Aedes vexans nipponii collected at a densely populated district of Seoul, and many insect-specific viruses, belonging to Baculoviridae, Nimaviridae, and Iflaviridae, are found in Culex quinquefasciatus and Culex tritaeniorhynchus gathered in China and Kenya (11,19). In this study, we identified six novel RNA viruses from three species of hymenopteran Pteromalidae parasitoids of the rice weevil S. oryzae. We characterized organization of  Table S2. Purple circles on the branches of each tree relate to 50 to 100% bootstrap support calculated from 1,000 replicates. viral genomes, their phylogenetic locations, and structures of viral nonstructural protein. These results indicate that the positive-strand and negative-strand RNA viruses could be simultaneously present in these three wasps.
Picornaviruses identified here are WWPSRV-1, WWPSRV-2, and AcPSRV-1. The organization of these viral genomes varies greatly. WWPSRV-1 possesses four ORFs in the genome, whereas WWPSRV-2 and AcPSRV-1 possess monocistronic and dicistronic genomes, respectively. In addition, these three viruses have the similar fluctuating distribution of sequencing reads, with the higher coverage area near 39 ends of the viral genome ( Fig. 2A). The 39 bias of the viral genome coverage is likely to be attributed to the expression of subgenomic RNAs (sgRNAs) (20). However, it also cannot be ignored that RNA degradation or shearing during sample handling and RNA extraction may also lead to a controllable 39 bias in short-read RNA-seq data (21). The analysis of PCR results and transcriptome data (Fig. 1A and B and Fig. 2A) confirms that WWPSRV-1 and WWPSRV-2 can be detected in these three wasps. The observations demonstrate that these two viruses may have the possibility of interspecies virus transmission. Just like Drosophila C virus (DCV), they can infect a much broader range of Drosophila species through horizontal transmission (22). Significantly, WWPSRV-2, which belong to family Iflaviridae of which some viruses can result in lethal infections in silkworms (23) and honeybees (24), can also be detected in the wasp host rice weevil. Moreover, it has been reported that an iflavirus Dinocampus coccinellae paralysis virus, found in Dinocampus coccinellae wasps which parasitize the lady beetle Coleomegilla maculata, can replicate in lady beetle cerebral ganglia and thus induce changes in lady beetle behavior such as tremors, gait disturbance, and limitations in movement (25). Accordingly, WWPSRV-2 may also play a role in the virus-wasp-weevil interaction. WWPSRV-1 shows high sequence similarity with Nasonia vitripennis virus, for which the incomplete genome was identified from cDNA libraries of parasitoid wasps Nasonia vitripennis (26). In the phylogenetic tree, WWPSRV-1 could cluster as a distinct lineage in Nora virus-related clade. In this clade, the most studied virus is Nora virus, which is a persistent nonpathogenic virus that infects various Drosophila species (27,28). However, Nora virus-infected D. melanogaster exhibits an increase in immune-related gene expression over time (29). AcPSRV-1 has a dicistronic genome with gene organization similar to that of dicistroviruses. Phylogenetic analyses of the RdRp domain revealed that, clustering with Himetobi P virus (detected in the small brown planthopper Laodelphax striatellus) as closest relative, AcPSRV-1 belongs to genus Triatovirus (30). In this genus, black queen cell virus, which also infects honeybees (hymenopteran insect), causes broad health problems with risk of colony sustainability (31).
The arrangement of the genomes of AcNSRV-1 and LdNSRV-1 follows the typical basic five-ORF pattern (39-N-P-M-G-L-59) of mononegaviral genomes (37). The mRNA of each gene can be synthesized directly from the viral genome (38). The high-coverage area of sequencing reads of AcNSRV-1 and LdNSRV-1 reveals that the N, P, M, and G may have a high expression level. Moreover, AcNSRV-2 contains an additional P6 protein between G and L proteins (Fig. 2B). The small additional proteins interposed between the typical basic protein were often identified in the plant rhabdoviruses (39). Some of these proteins' functions have been made clear, e.g., P3 acts as a movement protein in rice transitory yellowing virus (RTYV) and rice yellow stunt rhabdovirus (RYSV) and P6 acts as a systemic RNA silencing suppressor in RYSV (40)(41)(42). In addition, a highly conserved motif in the intergenic junction which could be grouped into TS, IS, and TI is predicted to be present in all three viruses. Moreover, the 39 and 59 end sequences of AcNSRV-1 are complementary and can form a putative panhandle structure thought to be involved in genome replication (17).
In conclusion, we have expanded the diversity of RNA viruses in the parasitoid wasps, and these RNA viruses may play an important role in the biocontrol of rice weevils. For better application of this information, much work should be accomplished. Just like PpNSRV-1 (12), it also should be confirmed whether these viruses can affect the parasitism rate and the life span of parasitoid wasps. Moreover, the virus transmission among wasps and their host weevil should also be investigated.

MATERIALS AND METHODS
Insect rearing. Rice weevils (S. oryzae) were initially collected from a cereal storehouse in Hangzhou, China in 2017. They were successively reared on wheat seeds in the laboratory at 30 6 1°C, 80% relative humidity, and 12-h light/12-h dark photoperiod. The mated adult beetles were provided with the fresh wheat seeds to lay eggs for a week (the seeds were replaced weekly). After 10 to 15 days, eggs developed into elder larvae used as hosts for the wasps. The colonies of A. calandrae, L. distinguendus, and T. elegans were initially collected from the parasitized rice weevil hosts of the same sites. The adults of the three wasp species were independently reared on 10% (vol/vol) honey at 25 6 1°C, 50% relative humidity, 12-h light/12-h dark photoperiod. For parasitism, the mated adults of the three wasp species were placed together with the wheat seeds that contained the host larvae. These parasitized beetle larvae were reared in the same conditions as were the wasp adults. The adults of A. calandrae or L. distinguendus emerged after 15 days. The adults of T. elegans emerged after 20 days. Each wasp colony was independently reared in the laboratory at least three generations before the experiments.
Next-generation sequencing, sequence assembly, and virus-related sequence discovery. Total RNA from 20 female adults of A. calandrae, L. distinguendus, or T. elegans was extracted using TRIzol (Invitrogen, CA, USA). The RNA quality was assessed using NanoDrop 2000 (ThermoFisher Scientific, Waltham, MA) and an Agilent 2100 bioanalyzer (ThermoFisher Scientific, CA, USA). RNA concentration was measured using an Invitrogen Qubit 2.0 Fluorometer (ThermoFisher Scientific, CA, USA). RNA samples (RNA integrity number . 8) were sent to Novogene company (Beijing, China) for the subsequent library construction. Sequencing libraries were generated using NEBNext Ultra RNA Library Prep kit for Illumina (NEB, USA) following the manufacturer's recommendations and index codes were added to attribute sequences to each sample. The clustering of the index-coded samples was performed on a cBot Cluster Generation System using TruSeq PE cluster kit v3-cBot-HS (Illumia) according to the manufacturer's instructions. After cluster generation, the library preparations were sequenced on an Illumina Hiseq platform and paired-end reads were generated. The low-quality, adaptor-polluted, and high content of unknown base (N) reads were removed. The clean reads were performed with de novo assembly using Trinity program (43). After assembly, the remaining contigs were annotated with BLASTX (E values # 1e25) in the NCBI nonredundant protein sequence databases to identify the virus-related sequences.
Viral genome sequencing and detection. Primers for viral genome sequencing and detection were designed based on the virus genome-like contigs (Table S1). Total RNA from adult female wasps or adult weevils was extracted using TRIzol. Single-stranded cDNA was synthesized from the RNA using the TransScript one-step gDNA removal and cDNA synthesis SuperMix kit (TransGen Biotech, Beijing, China). cDNA was used as a template for PCR. The terminal sequences of the viral genome were confirmed by 59 or 39 RACE according to the instructions of the SMART RACE cDNA amplification kit (Clontech, CA, USA). All amplified PCR products were cloned into pGEM-T Easy vectors (Promega, Beijing Branch, China) and sequenced.
To identify putative conserved transcription termination and initiation sequences of (2)ssRNA viruses (negative-sense single-stranded RNA), noncoding viral genome regions were analyzed using the open-source MEME suite 5.3.0 for the motifs analysis (48).
To identify the viral genomic RNA distribution in these three wasps' transcriptomes, clean RNA-seq reads from all three wasps were mapped to the complete genome of the six novel RNA viruses using the CLC Genomics Workbench 12.

SUPPLEMENTAL MATERIAL
Supplemental material is available online only. TEXT S1, TXT file, 0.1 MB.