Expression profiles and phylogenetic properties of venom gland‐specific viruses in some aculeate bees and wasps

To identify viruses and compare their abundance levels in the venom glands of hymenopteran species, we conducted venom gland‐specific transcriptome assemblies and analyses of 22 aculeate bees and wasps and identified the RNA genomes of picornaviruses. Additionally, we investigated the expression patterns of viruses in the venom glands over time following capture. Honeybee‐infecting viruses, including the black queen cell virus (BQCV), the deformed wing virus (DWV) and the Israeli acute paralysis virus (IAPV), were highly expressed in the venom glands of Apis mellifera and social wasps. This finding suggests that the venom of bees and wasps is likely to contain these viruses, which can be transmitted horizontally between species through stinger use. Apis mellifera exhibited an increasing pattern of abundance levels for BQCV, DWV, IAPV and Triatovirus, whereas the social wasp Vespa crabro showed increasing abundance levels of IAPV and Triatovirus over different capture periods. This suggests that the venom glands of honeybees and wasps may provide suitable conditions for active viral replication and may be an organ for virus accumulation and transmission. Some viral sequences clearly reflected the phylogeny of aculeate species, implying host‐specific virus evolution. On the other hand, other viruses exhibited unique evolutionary patterns of phylogeny, possibly caused by specific ecological interactions. Our study provides insights into the composition and evolutionary properties of viral genes in the venom glands of certain aculeate bees and wasps, as well as the potential horizontal transmission of these viruses among bee and wasp species.


Introduction
The characteristic feature of viruses is their omnipresence in all life forms, which makes them a potential threat to the health of organisms through their transmission (Yanez et al. 2020).Pathogen transmission to new hosts is a fundamental process in disease ecology.Viral transmission can be broadly categorized into two types: horizontal and vertical transmission (Cory 2015).Insects exhibit a vast diversity of viruses with different transmission strategies.Horizontal transmission refers to the transmission of infectious pathogens between individuals of the same generation.In insect-infecting viruses, horizontal transmission commonly occurs via the ingestion of contaminated food, sexual contact or vector-mediated transmission (Cory 2015).Various honeybee-infecting viruses (HVs) are known to be horizontally transmitted through food-borne transmission, topical or body contact with contaminated bees, venereal transmission and vector-mediated transmission (Yanez et al. 2020).Vertical transmission, on the other hand, involves the transmission of viruses to the next generation.In the case of honeybees, some viruses are vertically transmitted from queens or drones to their offspring through venereal infections or sperm (Yanez et al. 2020).
Nest parasitism and usurpation are common methods of nest foundation for eusocial hymenopteran species (Fisher 1993;Litte 1979;Matthews & Matthews 1979).Bumblebee and wasp queens utilize their stingers to fight over nest ownership, and the use of stingers for nest usurpation has been observed in orchid bees (Bohart 1970), sweat bees (Knerer & Atwood 1966), anthophorid bees (Eickwort & Abrams 1980), paper wasps (Turillazzi et al. 1990) and yellow jacket wasps (Reed & Akre 1983).Stingers are also used for fighting among recently emerged honeybee queens (Fisher 1993).Vespid wasp workers often engage in fights to establish dominance over other workers using their stingers (Spradbery 1973), indicating the importance of stingers in defense, nest ownership and reproduction for eusocial hymenopteran species.
Social wasps pose a problem as predators of honeybees, inflicting significant damage on honeybee colonies (Dejong 1979;Ifantidis 2003;Tzanakakis & Katsoyannos 1998;Wegner & Jordan 2005).Wasps obtain proteins from workers or larvae and carbohydrates from nectar or honey within honeybee hives (Bacandritsos et al. 2006).As a result, they consume bee broods, steal honey and destroy colonies (Mayer et al. 1987).Direct contact between wasps and honeybees, particularly through stinging or eating infected honeybees (Mazzei et al. 2018), can increase the likelihood of social wasps becoming infected by HVs through venom exchange.In contrast, bumblebees share close phylogenetic relationships with honeybees (Lockhart & Cameron 2001) and exhibit similar ecological, foraging and temporal behaviors (Ruiz-Gonzalez & Brown 2006).They also share floral resources, including nectar and pollen, which indirectly and horizontally transmit viruses to both species via food-borne transmission.Taken together, our hypothesis is that there might be an overlap of specific viruses between the viromes of venom glands in wasps and honeybees that have been shaped by venom interactions, but that the viromes of bumblebees might display distinct expression profiles or phylogenetic properties of viruses, compared with those of wasps and honeybees.
Most protein-coding messenger RNAs in eukaryotes contain a poly-A tail, which consists of a long chain of adenine nucleotides.The presence of polyadenylated RNAs allows for their selection via oligo-dT priming during reverse transcription, making the sequencing of polyadenylated RNAs one of the most common applications in RNA-seq (Hrdlickova et al. 2017).Additionally, the poly-A tail is recognized as a characteristic and functionally important feature of all picornavirus RNA genomes (Kempf & Barton 2015).Picornaviruses, such as black queen cell virus (BQCV), deformed wing virus (DWV), Israeli acute paralysis virus (IAPV), sacbrood virus (SBV) and Triatovirus are commonly detected in honeybees, whereas the presence of picornaviruses has not been widely reported in the venom gland of social wasps.Therefore, our focus has been on identifying and investigating these picornaviruses from the venom gland transcriptomes of certain aculeate bees and wasps.
In this study, we reanalyzed the venom gland-specific transcriptome data of several aculeate bees and wasps (Yoon et al. 2020) and performed a comparative analysis of their viromes using in silico methods.With the rearranged transcriptome data, we reassembled and annotated viral genes using de novo assembly and the Basic Local Alignment Search Tool (BLAST) search engine, estimated phylogenetic relationships between viral sequences and conducted a comparative analysis of the expression profiles of viruses isolated from the venom glands of various aculeate bees and wasps to validate our hypothesis.

Insect collection and total RNA extraction
Female social wasps, including Parapolybia indica, Polistes daehanicus, Vespa ducalis, Vespa mandarinia and Vespa velutina, were collected from Taehwa Mountain at Chugokri, Docheok-myeon, Gwangju-si, Gyeonggi Province, Republic of Korea, in September 2022.All collected hymenopteran species were immediately frozen in liquid nitrogen and stored in a deep freezer until further use.For each species collected from the same area, 10 venom glands were dissected using forceps (FONTAX, Prilly, Switzerland) and homogenized together in 200 μL of TRI reagent (Molecular Research Center, Cincinnati, OH, USA) using a glass tissue grinder (Duran Wheaton Kimble, Mainz, Germany) for 5 min at 25°C.The homogenates were transferred to 1.7 mL Axygen Maxy Clear Microtubes (Corning, New York, NY, USA) and vortexed for 15 s after the addition of 20 μL of bromochloropropane (Molecular Research Center).The mixture was incubated for 10 min at 25°C, followed by centrifugation at 12 000 × g for 15 min at 4°C.A total of 100 μL of the aqueous phase was transferred to new 1.7 mL microtubes (Corning), and the RNAs were then precipitated by adding 100 μL of isopropanol (Sigma-Aldrich, St. Louis, MO, USA) and stored at 25°C for 10 min.After centrifugation at 12 000 × g for 10 min at 4°C, the supernatant was removed and the RNA pellet was then washed twice by centrifugation at 7500 × g for 5 min at 4°C using 75% ethyl alcohol.Ethyl alcohol was removed, and the total RNA pellets were air-dried for up to 5 min and dissolved in 50 μL of nuclease-free water (NFW) (Thermo Fisher Scientific, Waltham, MA, USA).
Construction of metatranscriptomic in silico cDNA libraries and virome analysis in the venom glands of some aculeate bees and wasps RNA-seq short reads from the venom glands of 17 aculeate bees and wasps were obtained from previous studies (Yoon et al. 2020;Yoon et al. 2022).Furthermore, short reads obtained from newly extracted total RNA from venom glands of five social wasps, including P. indica, P. daehanicus, V. ducalis, V. mandarinia and V. velutina, were obtained using previously described protocols (Yoon et al. 2022).Low-quality short reads were removed using the fastp program (Chen et al. 2018), and quality-filtered short reads were used as inputs for the metaSPAdes program in the SPAdes assembler package (Prjibelski et al. 2020), with default parameters, to obtain contig sequences.Contigs shorter than 1 kb were removed to construct the metatranscriptomic in silico cDNA libraries.
The collected samples used in this study were specifically identified by comparing their mitochondrial cytochrome c oxidase subunit I (mtCOI) genes, obtained from their transcriptome data, against the mtCOI data available from the National Center for Biotechnology Information (NCBI) (Table S1).All of the samples exhibited identities ranging from 90% to 100%, except for Sceliphron deforme and a Sphecidae species, which were not found in the NCBI database.The high identities observed in the mtCOI genes may indicate a high level of confidence in species classification.
The raw data of transcriptome sequencing from the following species in previous studies and recently examined species have been uploaded to the NCBI Short Read Archive (Table S2).
The sequences of the metatranscriptomic in silico cDNA libraries were subjected to the discontiguous MEGABLAST program in the Blast2GO program (Conesa et al. 2005) using the Virosaurus database (release Apr 2020_4.1;Gleizes et al. 2020) to search for insect virus sequences.Sequences with a BLAST e-value lower than 1e À100 were further analyzed using BLASTN and the NCBI NT database (released April 2023).To quantify the abundance of metatranscripts, quality-filtered short reads from each sample were mapped to an in silico cDNA library using the Kallisto program (Bray et al. 2016) for calculating transcripts per million (TPM) values.The calculation of TPM is expressed as follows: TPM = (exon reads in gene/total exon reads) × 1 000 000.This formula is used to normalize gene expression levels, considering the length of the gene and the total number of reads in the sample.

Quantitative real-time polymerase chain reaction (qPCR)
Total RNA was extracted from the venom glands of two social wasps, V. crabro and Vl.flaviceps, and one bumblebee, B. ignitus, which were collected and immediately frozen in liquid nitrogen or frozen after being captured for 2 and 6 h at 25°C (Yoon et al. 2022).RNA extraction was performed using 200 μL of TRI reagent (Molecular Research Center) following the protocols described earlier.To eliminate DNA contamination, DNase I treatment (TaKaRa, Kyoto, Japan) was performed.Up to 20 μg of total RNA was mixed with 4 units of recombinant DNase I (TaKaRa), 10× DNase I buffer, 8 units of ribonuclease inhibitor (Promega, Madison, WI, USA) and up to 20 μL of diethylpyrocarbonate (DEPC)-treated water (Thermo Fisher Scientific).The mixture was incubated for 30 min at 37°C.Recombinant DNase I (TaKaRa) was inactivated by adding 1 μL of 0.5 mol/L ethylenediaminetetraacetic acid (EDTA) (Sigma-Aldrich) and incubating at 80°C for 2 min.
For RNA precipitation, a total of 5 μL of 3 mol/L sodium acetate (Sigma-Aldrich) and 250 μL of chilled ethyl alcohol were added to the DNase I-inactivated mixture, which was then stored at À80°C for 20 min.After centrifugation at 12 000 × g for 10 min at 4°C, the supernatant was removed and the pellets were washed with 70% chilled ethyl alcohol by centrifugation at 12 000 × g for 5 min at 4°C.The pellets were air-dried for up to 5 min and dissolved in 50 μL of nuclease-free water (NFW; Thermo Fisher Scientific).
The cDNA was synthesized using SuperScript IV Reverse Transcriptase (Thermo Fisher Scientific).Up to 5 μg of DNase I-treated total RNA was mixed with 0.5 μL of 50 μmol/L oligo-dT primer, 0.5 μL of 10 mmol/L dNTP mix (Thermo Fisher Scientific) and up to 6.5 μL of DEPC-treated water (Thermo Fisher Scientific).The mixture was incubated at 65°C for 5 min and then placed on ice for at least 1 min.The cDNA synthesis mix, composed of 2 μL of 5× SuperScript IV buffer (Thermo Fisher Scientific), 0.5 μL of 0.1 mol/L dithiothreitol (DTT) (Thermo Fisher Scientific), 0.5 μL of ribonuclease inhibitor (Promega) and 0.5 μL of SuperScript IV (Thermo Fisher Scientific) was gently mixed with the total RNA mixture and incubated at 50°C for 10 min.The reaction was terminated by incubation at 80°C for 10 min, and the samples were stored on ice.
As a result of the limited quantity of RNA available for confirming the copy numbers of various genes in the transcriptome data, three viruses, including DWV, IAPV and Triatovirus, along with a reference housekeeping gene (dimethyladenosine transferase) were chosen for analysis.Primers specific to viral genes were designed with comparable lengths, %GC content and amplicon sizes (Table S3).A total of 10 μL of PCR reaction mixture was composed of 10 ng of cDNA, 5 μL of TB Green Premix Ex Taq II (TaKaRa), 5 pmol/L of sense and antisense primers and 1 μL of distilled water.The relative transcription levels of the viral genes were calculated using the Pfaffl equation (Pfaffl 2001) and compared with the TPM values of the venom gland-specific transcriptome data.

Identification of viral contigs in the in silico cDNA libraries from venom gland metatranscriptomes
The BLAST search results identified six viruses, DWV, BQCV, IAPV, Himetobi P-like Triatovirus, SBV and an unclassified Iflavirus from the in silico cDNA libraries of venom glands.
BQCV has been detected in A. mellifera (GenBank accession no.OR416411), two solitary hunting wasps, S. deforme (GenBank accession no.OR416392) and Sphecidae sp.(GenBank accession no.OR416393), and one social wasp, V. analis (GenBank accession no.OR416394).The amino acid sequences of BQCV showed high sequence identities among the four hymenopteran species.The BQCV contig sequences encoded the 817 amino acids of BQCV ORF2 structural polyprotein, whereas the 50 region (corresponding to 1-430 aa) for S. deforme was not fully assembled (Figure S1).Other sequences from the four species (413-817 aa) were highly conserved.The BQCV sequences of S. deforme and Sphecidae sp. were identical to those of A. mellifera, whereas those of V. analis showed 99.8% sequence identity, with two amino acid mismatches at positions 4 and 352.The phylogenetic tree of the BQCVs matched the phylogeny of the mtCOI genes (Fig. 1A,B).
DWV contigs were detected in seven hymenopteran species: A. mellifera, Sphecidae sp., V. analis, V. crabro, V. mandarinia, V. velutina and Vl.flaviceps (Figure S2).However, only the DWV sequences from A. mellifera and Sphecidae sp.encoded a full-length DWV polyprotein of 2893 amino acids.The contigs from the other hymenopteran species showed partial assembly or deletions within the DWV coding sequences.The DWV sequences of the six species examined showed more than 98% sequence identity compared with those of A. mellifera.Among them, Sphecidae sp.exhibited the highest percentage of sequence identity (99.6%), whereas Vl. flaviceps showed the lowest percentage of sequence identity (98.7%), compared with that of A. mellifera.Phylogenetic analysis of the DWVs showed a similar pattern to that of the mtCOI genes, except for V. crabro and Vl.flaviceps (Fig. 1A,C).The positions of V. crabro and Vl.flaviceps in the DWV phylogeny were switched compared with their positions in the mtCOI phylogeny.This suggests that the DWV of each species possibly reflects the phylogenetic relationships among aculeate bees and wasps.
IAPV contigs were obtained from five hymenopteran species: A. mellifera, V. analis, V. crabro, V. velutina and Vl.flaviceps.These contigs were found to code 1950 amino acids of IAPV open reading frame 1 (ORF1) nonstructural polyprotein, except for Vl.flaviceps, which had a 50 unassembled region and a partial deletion (Figure S3).The IAPV from Vl. flaviceps showed the highest sequence identity (99.2%) with that of A. mellifera, whereas that of V. analis exhibited the lowest sequence identity (98.1%).Phylogenetic analysis of IAPV sequences clearly reflected the evolutionary pattern of the aculeate species, as it showed a matching pattern with the phylogenetic tree based on the mtCOI genes (Fig. 1A,D).
Himetobi P-like Triatovirus was found in the venom gland of five hymenopteran species.The Triatovirus contigs from A. mellifera, B. ignitus, Vl. flaviceps, V. mandarinia and V. velutina were not completely assembled sequences (Figure S4).These species showed partially assembled 50 regions of their Triatovirus sequences and exhibited highly conserved translated amino acid sequences within the 30 coding regions for the virus capsid protein.The Triatovirus sequences of three species, including Vl. flaviceps, V. mandarinia and V. velutina, clearly matched with those of A. mellifera, whereas only B. ignitus showed 97.6% sequence identity.In the phylogenetic analysis of Triatovirus, A. mellifera and B. ignitus did not cluster together as observed in the phylogeny of the mtCOI genes of aculeate species (Fig. 1A,E).Instead, B. ignitus was grouped with V. mandarinia in the Triatovirus phylogeny.The Himetobi P-like Triatovirus has been reported to be present in the muscles and guts of infected V. velutina thus far (Dalmon et al. 2019).However, it is likely that this virus represents an enzootic infection across numerous aculeate species.
An unclassified Iflavirus contig from the A. flavomarginatum venom gland showed high sequence identity with the diamondback moth Iflavirus (NC_03438).This contig contained the coding sequence of 2950 amino acids corresponding to the coding region of the diamondback moth Iflavirus and showed 98.1% sequence identity with the amino acid sequence of the viral polyprotein.As the larvae of lepidopteran species are often the prey for predacious hymenopteran species (Prezoto et al. 2019), the viruses might transmit to Vespidae species or vice versa.
The SBV contig was only found from the venom glands of Sphecidae sp.This contig was a 30 partial sequence of the virus and showed high sequence identity with the SBV of A. mellifera (AWK58136).This contig contained the coding sequence for 1374 amino acids of the SBV polyprotein, and the translated sequence showed a perfect match with the N-terminal region of the SBV polyprotein.SBV is a prevalent honeybee virus that is frequently detected in both brood and adult bees (Truong et al. 2023).However, it has only been identified in the venom glands of Sphecidae sp. and not in other species.This result could indicate the limited tissue tropism of SBV, but further research is necessary to validate this observation.B. consobrinus and B. ussurensis).These results suggest that the close ecological relationships between honeybees and wasps affect the formation of viral flora in some aculeate bees and wasps, which might contribute to the existence of HVs in these wasps.Among the examined species, BQCV exhibited relatively high abundance levels in A. mellifera and Sphecidae sp.(TPM value = 2992 and 898, respectively), indicating a broad host spectrum for BQCV (Fig. 2).Although relatively lower abundance levels of BQCV were observed in the venom glands of S. deforme and V. analis, the presence of BQCV in the four species suggests the possible passage of BQCV from honeybees to solitary hunting and social wasps and the possibility that wasps could serve as vectors and additional hosts for BQCV.Further studies are required to better understand the possibility of the vertical transmission of BQCV in solitary hunting and social wasps, as BQCV is vertically transmitted from the queen to her offsprings in honeybee colonies (Spurny et al. 2017).
DWV exhibited high abundance levels in the venom glands of A. mellifera (GenBank accession no.OR416412), four Vespa species, V. analis (GenBank accession no.OR416396), V. crabro (GenBank accession no.OR416397), V. mandarinia (GenBank accession no.OR416399) and V. velutina (GenBank accession no.OR416400), one Vespula species (Vl.flaviceps, GenBank accession no.OR416398) and one Sphecidae species (GenBank accession no.OR416395) (Fig. 2).Among these species, A. mellifera showed an extremely high abundance level of DWV (TPM value = 50 015), indicating that it is a potential reservoir for DWV.Additional studies should be conducted to elucidate the clinical evidence of DWV-infected wasps and the potential vertical transmission routes of DWV in solitary hunting and social wasps.IAPV was highly abundant in the venom glands of A. mellifera (GenBank accession no.OR416413) and three Vespa species, V. analis (GenBank accession no.OR416401), V. crabro (GenBank accession no.OR416402) and V. velutina (GenBank accession no.OR416404), and showed a relatively lower abundance level in the venom glands of Vl. flaviceps (GenBank accession no.OR416403) (Fig. 2).Among the examined species, V. analis and V. velutina exhibited extremely high abundance levels of IAPV (TPM values = 110 039 and 1126, respectively).The high abundance levels of the HVs, including BQCV, DWV and IAPV, in social wasps strongly suggest a close relationship between honeybees and wasps.
The Himetobi P-like Triatovirus was abundant in the venom glands of social wasps, including Vl. flaviceps (GenBank accession no.OR416406) and V. velutina (GenBank accession no.OR416408) (TPM values = 672 and 1423, respectively), but exhibited a relatively lower abundance level in Apidae species, including A. mellifera (GenBank accession no.OR416414) and B. ignitus (GenBank accession no.OR416405) (TPM values = 10 and 2, respectively) (Fig. 2).To investigate the sequence identities of Triatovirus, we used the BLAST search engine and observed that Triatovirus in the venom gland of five species, including A. mellifera, B. ignitus, Vl. flaviceps, V. mandarinia (GenBank accession no.OR416407) and V. velutina showed more than 81.6% sequence identity with the Triatovirus of A. mellifera, and less than 32.19% sequence identity with BQCV.These results indicate that the Triatovirus detected in these five species naturally infects hymenopteran species, and that the abundance pattern of Triatovirus exhibits close relationships between Apidae and Vespa species.
SBV was only observed in the solitary hunting wasp Sphecidae sp.(GenBank accession no.OR416409) with a low abundance level (TPM value = 4), whereas an unclassified Iflavirus was observed in A. flavomarginatum (GenBank accession no.OR416410) with a low abundance level (TPM value = 11) (Fig. 2).The partial sequence of SBV from Sphecidae sp.showed a clear match with that of the honeybee SBV (GenBank accession no.AWK58136.1).The Iflavirus from A. flavomarginatum showed a sequence similarity of 98.1% with the diamondback moth Iflavirus (GenBank accession no.QNS36256.1),and this is a case of an Iflavirus infecting both lepidopteran and hymenopteran insects, which is another example of putative viral transmission between prey and predator.
Hierarchical clustering analysis based on abundance levels revealed that HVs, including BQCV, DWV and IAPV, were clustered together, whereas Triatovirus, SBV and Iflavirus were clustered separately (Fig. 2).BQCV and DWV showed close relationships in the hierarchical tree structures.Apis mellifera exhibited a similar abundance pattern of viruses, including BQCV, DWV, IAPV and Triatovirus, to that of Vespa and Vespula species, such as V. crabro, V. mandarinia, V. analis, V. velutina and Vl.flaviceps.Notably, A. mellifera was clustered with the social wasp group, including V. analis and V. velutina, suggesting prospective virus exchanges between A. mellifera and Vespa species.

Expression patterns of viruses in the venom glands of some aculeate bees and wasps over different post-capture periods
In the previous study, the dynamic expression pattern changes of venom proteins in the venom glands of aculeate species under duress were demonstrated (Yoon et al. 2022), and the dynamics of viral RNA replications were therefore also assessed.The changes in abundance levels of four viruses in the venom glands of five hymenopteran species (A.mellifera, B. ignitus, P. rothneyi, V. crabro and Vl.flaviceps) exhibited six different patterns over different post-capture periods, as follows (Fig. 3): (i) decreasing and gradually increasing; (ii) constantly decreasing; (iii) decreasing and sharply increasing; (iv) increasing and gradually decreasing; (v) constantly increasing; and (vi) increasing and gradually increasing.The abundance level of the four viruses in the venom glands of A. mellifera has increased, depending on the time after capture, with the following three patterns: increasing and gradually increasing pattern in BQCV; increasing and gradually decreasing patterns in DWV and IAPV; and a constantly increasing pattern in Triatovirus (Fig. 3A).Abundance levels of DWV showed decreasing and gradually increasing patterns.IAPV and Triatovirus exhibited constantly increasing abundance levels in the venom glands of V. crabro over different post-capture periods (Fig. 3B).DWV exhibited decreasing and gradually increasing expression patterns, and IAPV and Triatovirus showed decreasing and sharply increasing expression patterns in the venom glands of Vl. flaviceps, depending on the different post-capture periods (Fig. 3C).Triatovirus exhibited constantly decreasing and decreasing and gradually increasing patterns in P. rothneyi and B. ignitus, respectively, depending on the post-capture periods (Fig. 3D,E).The higher average abundance level of HVs in A. mellifera and V. crabro after 6 h, compared with at the time of capture, indicates the possibility of viral replication in the venom glands.
DWV of V. crabro and Vl.flaviceps (Fig. 4A,D) and Triatovirus of B. ignitus (Fig. 4G) exhibited the highest TPM values at the time of direct capture, showing a decreasing expression pattern depending on the post-capture period.The relative transcription level of DWV from the two social wasps and Triatovirus from B. ignitus also showed a clear matching expression pattern with the TPM values (Fig. 4A,D,G).IAPV and Triatovirus of V. crabro (Fig. 4B,C) showed a gradually increasing expression pattern until 2 h after capture and a sharply increasing expression pattern at 6 h after capture.IAPV and Triatovirus of Vl. flaviceps (Fig. 4E,F) exhibited a decreasing expression pattern at 2 h post-capture and the highest abundance level at 6 h post-capture.The relative transcription levels of IAPV and Triatovirus of V. crabro and Vl.flaviceps also showed a clearly matched expression pattern with the TPM values (Fig. 4B,C,E,F).Similar expression patterns between the TPM values of the venom gland-specific metatranscriptome data and the relative transcription levels from the qPCR data imply that the estimated transcriptional abundance of viral genes from the metatranscriptome data is reliable.

Discussion
In this study, we identified various HV sequences via the de novo assembly of venom gland-specific transcriptome data and comparative analysis of viral contig sequences.As amino acid substitutions in viral proteins are known to affect their function and structural properties (Hanifa et al., 2023), we conducted a sequence identity analysis of amino acids in viruses using the BLASTX and TBLASTX programs.This analysis aimed to accurately determine the differences in viral sequences based on the host species of specific aculeate bees and wasps.The poor assembly of the 50 region or the absence of 50 partial sequences of HVs, such as those in the BQCV of S. deforme, IAPV of Vl. flaviceps and Triatovirus of A. mellifera, B. ignitus, Vl. flaviceps, V. mandarinia and V. velutina could be attributed to the presence of an internal ribosomal entry site (IRES) found in dicistroviruses.The picornaviral IRES, located upstream of the start codon, is known for its double-hairpin structure, which could pose challenges for sequencing (Hellen & Sarnow 2001).As the oligo-dT primers used for the reverse transcription of polyadenylated RNAs are known to generate high frequencies of truncated cDNA via internal poly-A priming (Nam et al. 2002), oligo-dT priming could also have contributed to the absence of 50 partial sequences of viral contigs.
BQCV and IAPV appeared to match the phylogenetic patterns of aculeate bees and wasps, whereas other viruses, including DWV and Triatovirus, did not.This suggests that not all viruses reflect the evolutionary pattern of aculeate lineages.These unique evolutionary patterns of viruses in aculeate bees and wasps may be influenced by the ecological traits of their host species, as suggested by the unique phylogenies of venom components in some aculeate bees and wasps (Yoon et al. 2020).In contrast to the hypothesis that adaptive venom traits provide selective advantages in certain hymenopteran species (Yoon et al. 2020), viruses may benefit from the ecological behavior of their insect hosts, leading to enhanced transmission efficiency and spread (Moreno-Delafuente et al. 2013).Therefore, viruses exhibited distinct evolutionary patterns in their phylogenetic trees.To validate this notion and obtain a clearer understanding of viral evolution, further studies on a larger scale are needed to explore the relationships between the phylogenetic properties of viruses and the ecological behavior of their insect hosts.Although HVs could replicate for several generations within a single host species, the entire sequence similarities of HVs in wasps demonstrated more than 97.6% identity compared with that of A. mellifera, indicating a high likelihood of inter-or intraspecies virus transmission in Hymenoptera.Additionally, the direct introduction of viral particles from venom fluid into the hemolymph of bees or wasps through their stingers could bypass the barrier of host specificities for viruses (Ekstrom & Hultmark 2016).
The venom gland of hymenopteran species, which is primarily composed of secretory epithelial cells, serves as a specialized organ for protein synthesis and secretion (Ferrarese et al. 2009).Consequently, this may provide viruses with a suitable environment for replication and secretion into the venom.Hymenopteran individuals of the same or different species could become infected by the venom containing secreted viruses during fights or self-defense.This suggests that virus transmission through venom could be one of the routes of inter-or intraspecies virus transmission in Hymenoptera.As viruses in the venom gland might be directly introduced into the hemolymph of bees or wasps through their stingers, this method represents the mechanical infiltration of the virus.Furthermore, it differs from the conventional infection route via the epithelial cells of digestive or reproductive organs.Consequently, viruses in the venom gland have the potential to bypass the tissue tropism barrier, which is important for primary infection and may lead to different outcomes of infection (Ekstrom & Hultmark 2016).
Analysis of the specific abundance of viruses in the venom gland is crucial for supporting the hypothesis that the venom glands of A. mellifera and various wasps might serve as reservoirs for viruses.This is because viruses could potentially reside within the cells of the venom gland, and given that melittin, one of the components of venom, has known antiviral properties, a high TPM value of viruses might not necessarily indicate the presence of infectious viral particles within the venom fluid.However, it is worth noting that arthropod-borne viruses (arboviruses) are recognized to replicate in the salivary glands of mosquitoes.This case serves as supporting evidence for the notion that the venom gland could indeed operate as a reservoir for viruses (Sanchez-Vargas et al. 2021).Therefore, conducting further investigations into the viral particle titers present in the venom fluids of Hymenopteran species is essential to provide additional support for this hypothesis.
Melittin is the main component of venom produced by the European honeybee A. mellifera, constituting 40%-60% of its dry weight (Liu et al. 2016;Memariani et al. 2020).This venom-derived peptide is composed of two connected α-helices (Terwilliger & Eisenberg 1982) and has been found to interact with cell membranes, induce pore formation at micromolar concentrations and cause cell lysis (Lee et al. 2013;van den Bogaart et al. 2008).Melittin possesses diverse biological activities, including antibacterial (Shi et al. 2016), antifungal (Do et al. 2014), anti-parasitic (Adade et al. 2013;Pereira et al. 2016), anticancer (Gajski & Garaj-Vrhovac 2013), anti-inflammatory (Lee & Bae 2016), antidiabetic (Hossen et al. 2017), anti-biofilm (Picoli et al. 2017) and antiviral (Uddin et al. 2016) properties.Melittin has been shown to exhibit antiviral activity against enveloped viruses, such as influenza A virus, vesicular stomatitis virus, respiratory syncytial virus and herpes simplex virus.It also has inhibitory effects on the replication of non-enveloped viruses, including enterovirus and coxsackievirus, which belong to the Picornaviridae family (Uddin et al. 2016).All members of the Picornaviridae family possess a positive-sense single-stranded RNA genome and a non-enveloped icosahedral capsid (Maclachlan & Dubovi 2010).During viral entry into host cells, small capsid proteins are released from the virus, which interact with cell membranes and facilitate the delivery of the viral RNA genome into the cytoplasm to initiate replication (Panjwani et al. 2014).In virus-cell attachment and virus-entry assays, it was observed that the virus titer in melittin-treated cells did not show significant differences compared with cells infected with the virus alone (Uddin et al. 2016).This suggests that melittin is unable to inhibit the initial steps of the virus life cycle, such as infection and replication, but instead inhibits virus infectivity by directly interacting with the virus prior to attachment to host cells.Therefore, if the concentration of melittin in the venom gland of A. mellifera is limited, non-inhibited viruses from venom-secreting cells in the venom gland may be able to infect other cells and proceed to replicate in the venom glands.These hypotheses might explain the exceptionally high abundance levels of HVs, such as BQCV and DWV, in the venom glands of A. mellifera.
Apis mellifera exhibited an increasing pattern of expression for all four viruses in the venom gland, whereas the social wasp V. crabro showed a constantly increasing abundance level of IAPV and Triatovirus over different post-capture periods.It is worth noting that after a 2 h capturing period of honeybees and wasps, venom extracts from the stinger were observed on the surfaces of the captured tubes.This observation could be attributed to the stressful environment for both honeybees and wasps, leading to the activation of venom glands and the production of venom extracts.The abundance level of viral fragments in the RNA-seq data may reflect the concentration of accumulated viral particles or replicated viruses in the venom gland.The observed increasing expression patterns over different post-capture periods of the viruses in these circumstances may suggest that the venom gland of honeybees and wasps provides a suitable environment for viral replication.Consequently, the RNA-seq data may reveal the abundance levels of replicated viruses.
The high abundance levels of HVs, including BQCV, DWV and IAPV, in the venom glands of A. mellifera, as well as in Vespa and Vespula species, such as V. analis, V. crabro, V. mandarinia, V. velutina and Vl.flaviceps, suggest the possibility of horizontal transmission between honeybees and social wasps.The transmission routes of these three HVs have been extensively studied.Horizontal transmission among adult honeybees (Singh et al. 2010) and vertical transmission from the queen to her offspring are considered the main reasons for the chronic persistence of BQCV in honeybee colonies (Spurny et al. 2017).DWV is known to be transmitted between honeybees through both horizontal transmission (fecal-oral) and vertical transmission (parent-offspring) (Chen et al. 2005;Chen et al. 2006;Yue & Genersch 2005).Infection with IAPV can be detected through vector-mediated and food-borne horizontal transmission pathways between honeybees, as well as vertical transmission from the queen to her progeny (Chen et al. 2014).Direct contact between honeybees and social wasps, particularly through stinging behaviors, may increase the likelihood of HV infection in the venom glands of Vespa and Vespula species.
Among the 22 aculeate bees and wasps examined, eight species, including A. mellifera, two solitary hunting wasps (S. deforme and Sphecidae sp.) and five social wasps (V.analis, V. crabro, V. mandarinia, V. velutina and Vl.flaviceps), were found to harbor HVs.One common characteristic of solitary hunting and social wasps is their predation on honeybees as a food source (Power et al. 2022;Strohm & Linsenmair 1997), suggesting that the ingestion of infected honeybees and larvae could serve as a possible horizontal transmission route for HVs (Mazzei et al. 2019).Although not all HVs were observed in all of the aculeate bees and wasps examined, the presence of HVs in the venom glands of solitary hunting and social wasps, and their absence in the venom glands of bumblebees, could be attributed to the direct contact facilitated by stinger use during predatory activities.Thus far, replicative forms of HVs, such as DWV (Mazzei et al. 2018), BQCV and the Kashmir bee virus (Mazzei et al. 2019), have been detected in V. velutina, and DWV has also been observed in V. crabro (Forzan et al. 2017).Most of these pathogens were identified through reverse transcription-polymerase chain reaction (RT-PCR) amplification using the entire body of individuals as a template.Positive-strand viruses such as HVs require the production of negative-strand RNA intermediates for replication (Yanez et al. 2012).Negative-strand RNA of IAPV has been detected in A. mellifera and in the head, thorax and abdomen of V. velutina, indicating viral replication in these two species (Yanez et al. 2012).The presence of replicative forms of HVs in social wasps and negative-strand RNA of IAPV in the abdomen of V. velutina supports the hypothesis that HVs could also replicate in the venom glands of social wasps.However, it is still unclear whether HVs can cause symptomatic infections or be transmitted vertically in social wasps.Therefore, further investigation of symptomatic infections, replication efficiencies and the transmission routes of HVs in social wasps is necessary to understand the role of HVs in invasive wasp species.

Conclusion
This study represents the first comparative systematic analysis of venom gland-specific metatranscriptome data, enabling the identification of expression profiles and phylogenies of viruses in 22 aculeate bees and wasps.Honeybee-infecting viruses, including BQCV, DWV and IAPV, were highly expressed in the venom glands of A. mellifera and social wasps, suggesting that the viruses could be transmitted horizontally between species through stinger use.An increasing pattern of abundance levels for HVs in A. mellifera and V. crabro over different capture periods indicates that the venom glands of honeybees and wasps may provide suitable conditions for active viral replication and may serve as an organ for virus accumulation and transmission.As a result, our study provides valuable insights into the composition and evolutionary patterns of venom gland-specific viral genes in some aculeate bees and wasps.

Figure 1
Figure 1 Phylogenetic analysis of mtCOI (A), black queen cell virus (B), deformed wing virus (C), Israeli acute paralysis virus (D) and Triatovirus (E) from the venom gland of some aculeate bees and wasps.

Figure 2
Figure 2 Hierarchical clustering of six viruses in 22 aculeate bees and wasps (left and top panel), and their differential expression profiles represented in heat map (right panel).Each row represents the viral genes and each column represents the species examined.Asterisks indicate that the honeybee-infecting viruses exhibited high expression levels in Apis mellifera and five social wasps, including Vespa analis, Vespa crabro, Vespa mandarinia, Vespa velutina and Vespula flaviceps.

Figure 3
Figure 3 Viral gene expression profiles of black queen cell virus, deformed wing virus, Israeli acute paralysis virus and Triatovirus in the venom glands of Apis mellifera (A), Vespa crabro (B), Vespula flaviceps (C), Polistes rothneyi (D) and Bombus ignitus (E) over different post-capture periods.The x-and y-axes represent the capture period and the log10 of the fold change in expression levels, respectively.The red solid lines in panels B, C and D indicate the mean values of gene expression levels.

Figure 4
Figure 4 Comparison of transcripts per million (TPM) values and the relative transcription levels of three viruses including deformed wing virus (A and D), Israeli acute paralysis virus (B and E) and Triatovirus (C, F, and G) from Vespa crabro, Vespula flaviceps and Bombus ignitus.

Table 1
Summary of the venom gland-specific metatranscriptome libraries in the venom glands of 22 aculeate bees and wasps Four HVs, including BQCV, DWV, IAPVand the Triatovirus, were highly abundant in A. mellifera.Three HVs were detected in the three social wasp species: V. analis (BQCV, DWV and IAPV), V. velutina (DWV, IAPV and Triatovirus) and Vl.flaviceps (DWV, IAPV and Triatovirus) (Fig.2).Two HVs were observed in one social wasp species, V. crabro (DWV and IAPV) and in one solitary hunting wasp species, Sphecidae sp.(BQCV and DWV).The solitary hunting wasp S. deforme and the bumblebee species B. ignitus possessed BQCV and Triatovirus, respectively.SBV and Iflavirus were only observed in two solitary hunting wasp species, Sphecidae sp. and A. flavomarginatum, respectively.None of the viruses were detected in the other eight social wasp species (P.indica, P. daehanicus, P. rothneyi, P. snelleni, P. varia, V. ducalis, V. dybowskii and V. simillima), one solitary hunting wasp species (E.decoratus) and three bumblebee species (B.ardens, †The largest contig size for which at least 25%, 50% and 75% of bases are contained in the assembled contigs.‡Theratio of bases composed of G or C.Venom gland viruses in bees and waspsEntomological Research 54 (2024) e12752