Metagenomic Analysis of Non-Apis Wild Insect Pollinators and Predators Shows Virus Threat to Their Survival

Background: Insect pollinators provide major pollination services for wild plants and crops, this is necessary for both agriculture and ecology, we should protect their population size and diversity as much as possible. Currently, we have known that the virus disease of honeybee can cause serious damage to bee colony, but we know little about the virus of other wild pollinating insects. Here we investigated the virus host wild insect pollinators as a reminder that they are facing virus threaten too. Methods: Transcriptome sequencing was used to investigate the viruses of Non-Apis insect pollinators and predators. Furthermore, host and replicability of several novel viruses that weakly related to honey bee viruses were determined by using reverse transcription-polymerase chain reaction (RT-PCR). Results: Three honey pathogenic virus, ve viruses host insect and twenty-six novel viruses were detected. Seven novel viruses showed weakly similarity to honeybee pathogenic viruses and ve of them were determined can be replicated their genomes in the corresponding host by detected complementary strands of viral genomes, suggests that they may be pathogenic to the corresponding host. Conclusion: So many novel viruses detected in wild insect pollinators and their predators indicates that we don’t have a clear understanding of the virus composition in pollinator insect and related species. The same living environment provides the conditions for virus cross infection, we should be alert to this situation in the protection of pollinating insects.


Background
Insect pollinators are necessary for most owering plants, and they play a key role in wild plant reproduction and food security [1,2]. Honey bee is the most popular pollinator because of its large population and pollination capability. However, other pollinators are also necessary to maintain ecological stability and to increase crop yield. Pollinator biodiversity is critical for pollination quality in agricultural productivity and conservation of the ecosystem [3,4]. The number of insect pollinators, including bees, is declining worldwide [5,6]. Several factors underlying pollinator decline can be mainly summarized as follows: heavy use of pesticides; worldwide spread of parasites, especially Varroa destructor mites; diseases caused by pathogenic viruses and other pathogens; monoculture cropping and plant biodiversity reduction; and competition between native and invasive species [5,7,8]. Adding further complexity to the issue, many of these factors act simultaneously on insect pollinators and can exert additive or even synergistic effects; for example, V. destructor is a vector of several honey bee viruses. High-quality food can help honey bee ght viruses. Some insecticides can increase the sensitivity of bees to viruses by affecting their immune mechanism [9][10][11][12][13][14][15].
Virus infection is a serious threat to honey bee, which can be infected by at least 25 types of viruses [16,17]. Of these, seven viruses responsible for severe diseases that threaten beekeeping are the deformed The application of next-generation sequencing technology in virus discovery has extended our understanding of viruses [19,20]. More than ten new honey bee pathogenic viruses have been detected by sequencing technology. These include Apis mellifera lamentous virus (AMFV), Lake Sinai virus (LSV), Varroa orthomyxovirus-1 (VOV-1), Apis mellifera bynyavirus-1 and − 2, Apis dicistrovirus, Apis mellifera avivirus, Apis mellifera nora virus-1, Apis mellifera rhabdovirus-1 and − 2, and Big Sioux River virus [17,21]. Furthermore, the discovery of new subtypes of known viruses, such as DWV-type B, DWV-type C, and several subtypes of LSVs, depends on sequencing technology [22][23][24]. Although we have accumulated much genome information on honey bee-related viruses, the evolutionary relationships of some honey bee viruses have not been expounded clearly. CBPV has two genome segments (segment 1: 3674 bp, segment 2: 2305 bp); its special genome structure and putative RNA-dependent RNA polymerase (RdRp) sequence have very low similarity to other viruses. One can only speculate that it may have an association with nodavirus and tombusvirus [25]. LSV was classi ed as Nodaviridae, but the LSVcontaining clade was distinct from the Nodaviridae family in an unrooted phylogenetic tree, and its several subtypes were detected in honey bee [22,26].
Virus research on insect pollinators has mainly focused on bees. In addition, compared with other factors threatening insect pollinators such as pesticides and food shortage, viruses are more easily ignored until before causing a disease [27]. Non-Apis insect pollinators are necessary for the ecosystem, but their virological investigation is often overlooked [3,28]. The honey bee pathogenic virus DWV is known to spread widely in wild insect pollinators, and other pathogenic viruses such as BQCV, SBV, and IAPV were detected in non-Apis hymenopteran species [29][30][31]. In this study, we used next-generation sequencing technology to investigate the viruses of non-Apis pollinators and their predators, hoping to expand our understanding of the viruses of insect pollinators and to provide some information for the protection of wild insect pollinator diversity.

Sample collection and processing
More than 50 species of pollinators and their predators in the eld were collected in Xiangshan, Beijing, from August to October 2020. All insects were collected on owers (insect pollinators) or next to owering plants (wasp-like insect predators). They were immediately stored in 5-mL centrifuge tubes in an ice box, and transferred in − 80 ℃ refrigerator until processing. Insects were individually ground into powder in liquid nitrogen and mixed into two pools for subsequent sequencing. Picture of all insect samples can be found in Supplementary Fig. 1 ~ Fig. 4.

Total RNA extraction and sequencing
All insect samples were ground separately in frozen condition and divided into two pools: pool A (pollinators) and pool B (predators). It is worth noting that we cannot determine the feeding habits of insects according to their morphology, and this grouping is preliminary. Sample composition of the pools can be found in Supplementary Table 1. Total RNA of the two mixed pools was extracted using the RNApure Total RNA Kit (RN0302; Aidlab Biotechnologies Co. Ltd., Beijing, China) according to the manufacturer's protocol. Sequencing libraries were generated using VAHTS mRNA-seq v2 Library Prep Kit for Illumina (NR601-01; Vazyme, Nanjing, China) following the manufacturer's recommendations, and index codes were added to attribute sequences to each sample. The libraries were sequenced on an Illumina NovaSeq platform to generate 150-bp paired-end reads according to the manufacturer's instructions. Finally, we obtained a total of 37 Gb clean reads, data is stored in the SRA database with accesssion: SRR14554108 -SRR14554111.

Sequence assembly and virus detection
Sequencing reads were de novo assembled using Trinity v2.12.0 [32]. The contigs were compared to the GenBank NR database by Diamond v2.0.8.146 [33]. We used a python script to lter out non-virus contigs. Virus-like contigs were mapped using Bowtie 2 v2.4.2 [34] to evaluate the contig quantity. Virus ORFs were annotated based on the result of ORF nder v0.4.3 [35] and the structure of the most closely related viral genome. Conserved domains of new viruses were identi ed using NCBI CDD BLAST v2.11.0 [36]. We ignored viruses for which the whole genome could not be obtained. Sequences of these viruses were either fragmentary genomes or had very few reads.

Phylogenetic analyses
RdRp regions identi ed by CDD were aligned between target viruses and related viruses using MUSCLE v3.8.31 [37,38]. All ambiguously aligned regions were subsequently removed using TrimAl v1.2rev59 [39], The best-t model of amino acid substitution in each dataset was determined using ModelTest [40].
Phylogenetic trees were then inferred using the maximum likelihood method implemented in RaxML with 1000 bootstrap replicates [41]. Phylogenetic trees were displayed and annotated using FigTree v1.4.4.

Identi cation of virus hosts
We extracted RNA from each insect powder using the RNApure Total RNA Kit (Aidlab). cDNA was prepared using TaKaRa PrimeScript™ RT Reagent Kit (Perfect Real Time RR037A; Takara Biomedical Technology Co. Ltd., Beijing, China). The presence of a virus was determined using PCR, which was carried out with the rst-strand cDNA products using the 2×TSINGKE Master Mix (blue) (TSE004; Beijing TsingKe Biotech Co., Led., Beijing, China) in 25-µL reactions and speci c primers. The speci c primers were designed based on the assembled viral genome sequences. PCR products were con rmed by agarose gel electrophoresis and Sanger DNA sequencing. Details of primer information are shown in Supplementary Table 1.

Prediction of potential pathogenicity of novel viruses
We used the reverse complementary strand of the virus to predict the pathogenicity of novel viruses. The detection process and reagents used were the same as those in reverse transcription-polymerase chain reaction (RT-PCR), except for speci c primers. Speci c primers were designed based on the assembled viral genome sequences, and a tag sequence was added to the forward primer. In reverse transcription, RNA was reverse transcribed with tagged forward primers and the control group with random primers.
Then the tag was used as the forward primer and normal reverse primer for PCR detection of cDNA products. PCR products were also con rmed by agarose gel electrophoresis and Sanger DNA sequencing. Details of primer information are shown in Supplementary Table 1.

Detection of known viruses
To obtain the virus information from the sequencing data, we used the Trinity software to assemble the sequencing data, and then the contigs were annotated using the Diamond program before querying the GenBank non-redundant (NR) database. We detected eight known viruses in the annotated results ( Table  1). Three of these were the honey bee pathogenic viruses DWV, ABPV, and CBPV. May eld virus 1 was rst detected in Bombus terrestris in Lebanon and UK, Vespa velutina associated acypi-like virus was rst detected in V. velutina nigrithorax in France, Scaldis River bee virus was rst detected in Osmia cornuta in Belgium, Arboretum almendravirus was rst detected in the mosquito Psorophora albigenu in Peru, and Hubei diptera virus 6 was rst detected in Diptera in China. The repeated detection of these viruses in different places and times indicates that they may be prevalent in their corresponding hosts and may exhibit pathogenicity.  (3). The two unclassi ed viruses were related to CBPV and Negevirus ( Fig. 1 and Table 2).

Identi cation of hosts of novel viruses
Because the pool was mixed with several species before sequencing, the hosts of novel viruses could not be obtained from the sequencing data. We selected seven viruses (i.e., XOLV, XTLV, XPLV2, XIV, XNLV, XPLV1, and XPLV4) to identify their hosts. These viruses were presumed to have a segmented genome or were weakly related to key honey bee viruses (i.e., LSV, CBPV, SBV, and DWV). Primers were designed based on novel viral genomes, and these primers were used to screen all samples before RT-PCR to con rm the sample source of novel viruses. Hosts of the novel viruses are listed in Table 2 and shown in Fig. 2. For segmented genome viruses, we further designed speci c primers for each segment to con rm that these segments come from one sample. All fragments showed positive results in the corresponding host samples (Fig. 3). Single-stranded RNA (ssRNA) viruses produce reverse complementary strands during replication. Whether an ssRNA virus can replicate in the host is determined by the presence of the reverse complementary strand of the ssRNA virus in the host [42][43][44][45]. So, to investigate if a novel virus replicated in the host, we screened for the presence of the positive-sense RNA strand of negative-sense genome viruses and the negative-sense RNA strand of positive-sense genome viruses using RNA-strand sense-speci c primertagged RT-PCR (Fig. 4). Reverse complementary strands of XOLV, XPLV1, XPLV2, XNLV, and XPLV4 were detected in the corresponding host samples. Speci c primer information can be found in Supplementary  Table 2.

Discussion
RNA sequencing (RNA-Seq) is a powerful tool that has been used to discover novel RNA viruses and to detect pathogenic viruses [19,46,47]. Many studies focus on the viromics of honey bee using sequencing methods, and this has greatly improved our understanding of bee viruses [16,22,[48][49][50]. Compared to honey bee, there are few virus studies on non-Apis insect pollinators, although their survival is under similar threat as honey bee [6, 31,51,52]. By using RNA-Seq technology, we preliminarily explored the viral composition of non-Apis pollinators and their predators.
Honey bee pathogenic viruses spread to non-Apis wild insect pollinators, and at least seven viruses (DWV, BQCV, SBV, IAPV, ABPV, SBPV, and CBPV) were detected in non-Apis insect pollinators. DWV and BQCV were able to replicate their genomes in bumble bees [29,31,45]. DWV, ABPV, and CBPV were detected in our sequencing data, and our samples did not include A. cerana, A. mellifera, or other Apis species, indicating that the three honey bee pathogenic viruses were detected in non-Apis wild insect pollinators.
Our results were consistent with previous reports. Spread of honey bee viruses to other pollinators should be noticed, especially when bees are used to provide commercial pollination services [53]. There is no doubt that commercial pollination services have created great value for agricultural production [54][55][56], but improving virus detection is an effective means to prevent the transmission of viruses from commercial insect pollinators to wild insect pollinators.
Five previously reported insect viruses (host in parentheses) were found in our samples. May eld virus 1 (B. terrestris), V. velutina associated acypi-like virus (V. velutina nigrithorax), Scaldis River bee virus (O. cornuta), Hubei diptera virus 6 (Diptera), and Arboretum almendravirus (P. albigenu) have been previously reported in other places and times [20,[57][58][59], and these viruses were also detected in our samples, indicating that they may be prevalent in the corresponding host and may exhibit pathogenicity.
We detected 26 novel RNA viruses by assembling sequencing data. Bee-infecting viruses are primarily positive-sense ssRNA (+ ssRNA) viruses of the order Picornavirales [27,60]. Seven of these novel viruses were classi ed in the order Picornavirales, suggesting that viruses of Picornavirales also infect non-Apis insect pollinators. Viruses in the families Rhabdoviridae, Flaviviridae, Orthomyxoviridae, and Sinhaliviridae can infect honey bee [16,21]. Novel viruses belonging to these families were detected in this study, suggesting that non-Apis insect pollinators and their predators have similar virus classes as honey bee. But the difference is that ve rhabdoviruses were detected in this study and only two rhabdoviruses have been reported in honey bee. This clearly infers that rhabdoviruses are more common in non-Apis pollinators and predators. Virgaviridae consists of plant viruses [61], and three novel viruses from this family were detected in this study. We think that these three viruses may have come from the contamination of plants. Due to the complex living environment of insect pollinators, we cannot exclude novel plant, fungal, and bacterial viruses. Seven novel viruses (i.e., XPLV2, XPLV4, XPLV5, XPLV1, XNLV, XIV, and XOLV) showed amino acid similarity to honey bee pathogenic viruses (i.e., DWV, SBPV, SBV, LSV, CBPV, and VOV-1). Five of these could replicate their genomes in the corresponding host, suggesting that they are pathogenic to the corresponding host. In addition, the hosts of these viruses closely interact with honey bee in the same ecosystem, suggesting that they are a potential threat to bees, just as host honey bees may also spread the virus to non-Apis insect pollinators [30].
The decline in insect pollinators is simultaneously driven by several factors and could act synergistically [5,62]. The interaction of multiple factors makes it more di cult to understand the reasons for the decline in insect pollinators. Some insecticides can exert additive or synergistic effects on virus-induced mortality and replication in honey bees by affecting their immune system [12]. Diet and virus infection are related in honey bee. High-quality diet has the potential to reduce mortality in the face of infection with IAPV [10]. Base on this study, we can see that the threat of the virus to wild pollinator insects may exceed our expectations. Not only the spread of bee virus threatens them, but also the virus carried by them may pose a greater threaten. virus research on insect pollinators mainly focuses on bees, but research on non-Apis virus host species and have same habitat as them is scarce.

Conclusions
In conclusion, the study of insect pollinators virology should expand the scope of species, there is a large gap in viruses in wild pollinators and predators that live in the same environment as honey bees, and corss host transmission of the virus should be noted in the future.

Declarations
Ethics approval and consent to participate Not applicable.

Consent for publication
The authors declare that they agreed to publish this paper with the permission of the publishing houses.

Availability of data and material
This ~550 bp sequence of PBNSPaV from grapevine in this study was submitted to GenBank with the accession number MH371356.

Competing interests
The authors declare that they have no competing interests. Authors' contributions WL conceived and designed the study. YH and DW performed the experiments. NL analyzed the epidemiological data. NL wrote the manuscript. All authors read and approved the nal manuscript. Figure 1 Phylogenetic relationship of novel viruses inferred from conserved RdRp amino acid sequences. RdRp conserved sequences of novel viruses and related viruses in the trees were identi ed by CDD. Trees were constructed using the maximum likelihood method, and the amino acid substitution model is annotated below each tree. The orthomyxo-like virus tree was based on polymerase subunit PA amino acid sequences, and the narna-like virus tree was based on complete RdRp amino acid sequences.

Figure 2
Hosts of seven novel viruses. Sample number is shown below each host.