Expanding diversity of bunyaviruses identified in mosquitoes

Mosquitoes interact with various organisms in the environment, and female mosquitoes in particular serve as vectors that directly transmit a number of microorganisms to humans and animals by blood-sucking. Comprehensive analysis of mosquito-borne viruses has led to the understanding of the existence of diverse viral species and to the identification of zoonotic arboviruses responsible for significant outbreaks and epidemics. In the present study on mosquito-borne bunyaviruses we employed a broad-spectrum RT-PCR approach and identified eighteen different additional species in the Phenuiviridae family and also a number of related but unclassified bunyaviruses in mosquitoes collected in Zambia. The entire RNA genome segments of the newly identified viruses were further analyzed by RNA sequencing with a ribonuclease R (RNase R) treatment to reduce host-derived RNAs and enrich viral RNAs, taking advantage of the dsRNA panhandle structure of the bunyavirus genome. All three or four genome segments were identified in eight bunyavirus species. Furthermore, L segments of three different novel viruses related to the Leishbunyaviridae were found in mosquitoes together with genes from the suspected host, the Crithidia parasite. In summary, our virus detection approach using a combination of broad-spectrum RT-PCR and RNA sequencing analysis with a simple virus enrichment method allowed the discovery of novel bunyaviruses. The diversity of bunyaviruses is still expanding and studies on this will allow a better understanding of the ecology of hematophagous mosquitoes.


Detection of diverse bunyaviruses within the Phenuiviridae and unclassified bunyaviruses
Field mosquitoes were collected from various locations in Zambia from 2014 to 2022.A total of 17,281 adult female mosquitoes were divided into 964 pools according to species for the purpose of detection and isolation of bunyaviruses (Supplemental Table S1).The broad-spectrum RT-PCR assay designed to detect a wide range of arbovirus species within Phenuiviridae was applied to RNA extracted from the mosquitoes and this allowed the detection of bunyavirus sequences from 150 pools 9 .BLAST analyses of sequences of the PCR products identified various bunyavirus-like sequences with some similarity to previously-identified viral RNA-dependent RNA polymerase (RdRP) genes (Table 1).

Regional difference of Cx bunyavirus 2 prevalence in Culex quinquefasciatus in Zambia
Among the 150 bunyavirus-positive pools, sequences of the PCR products from 103 Culex quinquefasciatus (Cx.quinquefasciatus) pools were closely related to Cx bunyavirus 2 (Accession No. MH188052) in the unclassified bunyavirus group with 87% nucleotide identity in the partial RdRP gene.As these new strains of Cx bunyavirus 2 were highly prevalent in Cx. quinquefasciatus, we further analyzed regional differences of the prevalence of Cx bunyavirus 2 infection (Table 2 and Fig. 2).The Cx bunyavirus 2 positivity was especially high in Cx. quinquefasciatus present in Livingstone (79/85 pools, minimum infection rate per 1000: 33.64-40.89)and Isoka (7/9 pools, minimum infection rate per 1000: 75.27) districts compared to other regions.The Cx. quinquefasciatus mosquito is the most common mosquito species in Zambia and indeed high numbers of Cx. quinquefasciatus were collected in all regions of the country.Although environmental factors such as rainfall, temperature, and animal populations that are characteristic of Livingston or Isoka have not been identified 9 , the obtained results suggest that regionally different factors may be involved in the transmission of Cx bunyavirus 2.

Viral RNA enrichment for mosquito RNA sequencing by RNase R treatment
Virus isolation assays of these viruses were attempted using several mosquito and mammalian cell lines, however, no productive viral growth was observed.To identify the entire tripartite RNA genomes of the detected bunyaviruses, we performed RNA sequencing analysis of the total RNA from the virus-positive mosquito lysates with a viral RNA enrichment step (Fig. 3A).In addition to the digestion of DNAs with DNase I, total RNA was treated with RNase R, which is a 3' to 5' exoribonuclease that digests linear RNAs, including cellular mRNAs, but does not digest lariat or circular RNA structures or double-stranded (ds) RNA with 3' overhangs shorter than seven nucleotides.Therefore, the bunyavirus genome, in which the 3' and 5' ends of each genome segment are highly complementary and form a dsRNA stem structure 15 , are expected to be resistant to RNase R treatment.First, we compared the number of bunyavirus reads from five representative RNA samples treated with or without RNase R.Although three segments of bunyavirus contigs longer than 1000 bases were obtained from both DNase I alone and DNase I and RNase R-treated samples, the number of reads on the bunyavirus contigs increased in the samples treated with RNase R compared to that with DNase I alone (Fig. 3B, C).These results indicate that RNase R treatment can enrich bunyavirus genes from total RNAs consisting of mostly host mosquito RNAs and this increases the reads of viral gene sequences (Supplementary Fig. S1).www.nature.com/scientificreports/

Identification of bunyavirus genome segments with conserved protein coding sequences
The RNA sequencing and BLASTx analyses eventually identified all three or four segments of the virus genome in eight bunyavirus species.As a result of attempting to identify the terminal sequences of these bunyaviruses by RACE analysis, we were able to determine the complete genome sequences only for Cx tenuivirus, Cx hudovirus and Cx bunyavirus 2 (Fig. 3D).The phylogenetic analyses of the putative coding sequences of RdRP on the L segment, glycoprotein on the M segment, and nucleoprotein sequences on the S segment were performed (Fig. 4A-C).RdRP and nucleoprotein sequences of Cx tenuivirus clustered with sequences of the genus tenuivirus, which has four segments of viral genome and infects host plants and vector insects.The putative glycoprotein on the M segment was found in Cx tenuivirus as well as previously identified Fitzroy Crossing tenui-like virus 1 16 , although glycoproteins of previously detected tenuiviruses have not been identified.The phylogenetic analysis of bunyavirus glycoproteins showed that Culex tenuivirus and Fitzroy Crossing tenui-like virus 1 have unique glycoprotein-like proteins distinct from the other bunyaviruses (Fig. 4B).In some samples, only the L and/or S segment of bunyavirus were identified by BLAST analysis.Three different L segments of bunyaviruses related to the Leishbunyaviridae (Cx leishbunyavirus 1, 2, and 3) and a putative S segment (Cx phenuivirus 3) were discovered in the Culex mosquito RNA samples.As bunyaviruses in the Leishbunyavirudae have been identified from protozoa in the family Trypanosomatidae 17 , we further explored protozoa genes in the RNA sequencing data in which Cx leishbunyavirus genes were found.As expected, contigs of Crithidia fasciculata ribosomal DNA gene were detected in the pool of Cx. quinquefasciatus (2015_Livingston_#15, Table 1) that contains two species of leishbunyaviruses, suggesting that these leishbunyaviruses may be derived from Crithidia fasciculata (Supplementary Table S2 and S3).Additionally, in the bunyavirales, RNA sequencing analyses in this study identified the L segment of previously known unclassified bunyaviruses, Kiristianstad virus from a pool of Cx. quinquefasciatus (18mong_50), and all three segments of Herbert virus in the Peribunyaviridae from a pool of Cx. quinquefasciatus (16mwi_21) (Fig. 4A-C).

Discussion
This study has revealed the existence of diverse bunyaviruses which have adapted to each mosquito species.The broad-spectrum RT-PCR assay designed to target the highly conserved region of RdRP gene within phleboviruses, including RVFV, allowed the detection of a wide variety of phenuiviridae-related viruses, including those plant-associated viruses.Thus, the broad-spectrum RT-PCR targeting the conserved RdRP gene is convenient for comprehensive detection of known arboviruses and unknown viruses in vector mosquitoes, ticks 18 , and other biological samples.The transmission cycle of the identified viruses is still unknown.Further investigation of animals and environmental samples and isolation of the viruses are required to clarify the individual viral life cycles, whether mosquitoes are a reservoir or just a transient vector, and the zoonotic potential of newly identified bunyaviruses.One of the discovered viruses, Cx tenuivirus with characteristic four genome segments, is closely related to the Genus tenuivirus, suggesting a possible plant host.When the adult mosquitoes feed on plant nectars, viruses can be transmitted horizontally between mosquitoes and plants.Interestingly, we found that Cx bunyavirus 2 was highly prevalent in Cx. quinquefasciatus which inhabit Livingstone and Isoka districts compared to other regions in Zambia.It is possible that the presence of other hosts of Cx bunyavirus 2, which is involved in the life cycle of Cx. quinquefasciatus, may be associated with the regional differences in virus positivity.The Livingstone City has a tourist spot of The Victoria Falls World Heritage Site/ Mosi-oa-Tunya National Park.A notable mosquito ecology in Livingstone which is the main source of Cx. quinquefasciatus appears to be a large waste stabilization pond covered with exotic water hyacinths adjacent to the National Park.It is likely that Livingstone has an ecological niche that contain non-mosquito hosts for Cx bunyavirus 2, and possibly other aquatic organisms in the aquatic stage of mosquito larvae.www.nature.com/scientificreports/In total RNA sequencing analysis, it is often difficult to detect very small amounts of viral genes in host total RNAs.The poly-A selection method is not applicable for enrichment of viral RNAs except for viruses with a poly-A tailed genome.Commercially available ribosomal RNA (rRNA) depletion systems target human, mouse, rat, and bacteria, and therefore species-specific probe preparations are required to deplete the rRNA from the total RNA of non-model animals.The results of RNA sequencing showed that our RNase R treatment system is very useful for the enrichment of bunyavirus genomes which have a dsRNA panhandle structure at the 3'-and 5'-ends of each genome segment, although it is difficult to rule out the possibility that the RNase R resistance may involve not only the double-stranded structures but also viral RNA protection by viral nucleoproteins or binding of the L protein to the 3' and 5' termini 19,20 .In addition, we found that the RNase R treatment system enriched  The maximum likelihood-based phylogenetic tree was constructed using partial RdRP amino acid sequences of representative bunyaviruses and identified genes in this study.Bootstrap values higher than 50 are shown adjacent to the tree branches.The tree is drawn to scale with branch lengths representing the number of substitutions per site.Black circles indicate bunyaviruses identified by the broad-spectrum RT-PCR.Bunyavirus contigs identified from RNA sequencing analyses are marked with red circles.The organisms in which the viruses were identified are color-coded with a caption in the box at the right lower corner of the Figures.The viruses considered to be mammalian pathogens or plant pathogens are highlighted in light orange and light green colors, respectively.not only bunyaviruses but also flaviviruses, which are positive sense ssRNA virus with higher order structures at their RNA 3'-ends 21 , and dsRNA viruses such as reoviruses.However, reads of ssRNA viruses with poly-A tail in their genome, such as ifraviridae, were decreased by RNase R treatment as expected (Supplementary Fig. S1).These results indicated that viral RNA genomes which have double-stranded or higher order structures at their RNA 3'-ends are enriched by RNase R digestion.Thus, the RNase R treatment approach has the advantage for the detection of RNA viruses, such as bunyavirus, flavivirus, and reovirus family members, in a variety of samples.It is expected to further improve the efficiency of virus detection by combinations of this simple and low-cost method with other virus enrichment methods and virome analysis pipelines to efficiently detect unknown viral genomes 10,[22][23][24] .For some bunyaviruses, identification of all genome segments was unsuccessful.Besides a low copy number of viral genome in mosquitoes, it is possible that the viral genome could not be recognized by BLAST analysis due to low homology with known viral genes whose sequences have been registered.In particular, the M and S segments are less conserved than the RdRP-coding L segment, making them difficult to identify.For instance, the predicted nucleoprotein sequence (266 aa) of the Cx phenuiviurs 3, whose L and M segment were not found, have only 35% amino acid identity with previously known bunyaviruses.Although tick-borne phleboviruses without a M segment have been reported 25 , virus isolation is necessary to identify the characteristics of such novel bunyaviruses.The M and S segments of Cx leishbunyavirus 1, 2, and 3 were also missing, while previously reported leishbunyaviruses from Crithidia spp.have putative M and S segments 17 .In order to elucidate the relationship between these leishbunya-like viruses, Crithidia parasite, and Culex mosquitoes, and their role or pathogenicity, isolation of virus together with Crithida parasite will be absolutely necessary.Some phenuiviruses have an ambisense coding strategy and encode nucleoproteins and nonstructural proteins (NSs) on the S segment.Short open reading frames (approximately 100 aa) in the antigenomic direction of the S segment genes were also found in some of the identified bunyaviruses.However, amino acids sequences of these open reading frames have no homology with the NSs of phenuiviruses, and BLAST analysis failed to find any similar or related proteins.
The newly discovered viruses identified in this study indicate that field-collected mosquitoes harbor highly divergent bunyaviruses..We established three to five CDC light traps (John W. Hock Co., USA) with yeast fermentation to supply CO 2 , and three BG-sentinel traps (Biogents AG, Germany) at different locations in inhabited area.The traps were set in the afternoon and left until the following morning, over a period of five nights on average in each district.After species identification, one to up to 40 female mosquitoes were pooled from each species and stored at −80 °C.Mosquitoes were first identified morphologically with reference to the African mosquitoes identification keys [27][28][29] .For mosquitoes that could not be identified morphologically, DNA extracted from the mosquitoes was genetically confirmed via sequencing of the cytochrome oxidase I (COI) gene 30 of mosquito DNA as a standard DNA barcoding for molecular identification.

Screening of phenuivirus genes by broad-spectrum RT-PCR
Pooled mosquitoes were homogenized in Minimum Essential Medium containing 2% fetal bovine serum using the BioMasher (Nippi, Japan).RNA was extracted from 100 µL of supernatants of mosquito homogenates using the Direct-Zol kit (Zymo research, USA) according to the manufacturer's instructions.The remaining mosquito homogenates were filtered and inoculated onto cells for virus isolation.To detect multiple phleboviruses, we designed degenerate primer sets targeting the conserved RdRP gene region; L-2779F 5'-CAR CAT GGW GGT YTDAGR GAR ATCTA-3' and L-3287R 5'-TGCARKATKCCY TGC ATCATHCCWG-3′ 9 .RNA samples were amplified using a PrimeScript One-step RT-PCR kit Ver.2 (Takara, Japan) and 1 µM of primer sets.The cycling protocol was comprised of 30 min of incubation at 50 °C for cDNA synthesis, followed by 2 min of incubation at 94 °C, 43 cycles each of 94 °C for 30 s, 52 °C for 30 s and 72 °C for 30 s, and 72 °C for 5 min.The PCR products were sequenced using a BigDye Terminator v3.0 Cycle Sequencing kit on an ABI PRISM 3130 Genetic Analyzer (Applied Biosystems, USA).

Molecular phylogenetic analysis
Deduced amino acid sequences of pan-phlebovirus RT-PCR products, or coding regions of the L, M or S segments were aligned with previously-characterized bunyaviruses (Supplementary Table S4) using the ClustalW 31 .

Figure 1 .
Figure1.Molecular phylogenetic analysis of partial RdRP sequences.The maximum likelihood-based phylogenetic tree was constructed using partial RdRP amino acid sequences of representative bunyaviruses and identified genes in this study.Bootstrap values higher than 50 are shown adjacent to the tree branches.The tree is drawn to scale with branch lengths representing the number of substitutions per site.Black circles indicate bunyaviruses identified by the broad-spectrum RT-PCR.Bunyavirus contigs identified from RNA sequencing analyses are marked with red circles.The organisms in which the viruses were identified are color-coded with a caption in the box at the right lower corner of the Figures.The viruses considered to be mammalian pathogens or plant pathogens are highlighted in light orange and light green colors, respectively.
www.nature.com/scientificreports/Methods Mosquito collection Mosquito collections were conducted between 2014 and 2022 with permission from the Department of National Parks and Wildlife, Ministry of Tourism and Arts of the Republic of Zambia, Excellence in Research Ethics and Science (ERES) converge ethics committee (IRB No: 00005948), and University of Zambia Biomedical research ethics committee (REF.NO. 1382-2020)

Figure 2 .
Figure 2. Prevalence of the Culex bunyavirus 2 in the Culex quinquefasciatus at each district in Zambia.The pie charts indicate the number of Culex bunyavirus 2-positive pools (yellow)/the pool number of Culex quinquefasciatus collected in each district shown on the map of Zambia.

Figure 4 .
Figure 4. Molecular phylogenetic analysis of the coding sequence of viral proteins.Maximum likelihood-based phylogenetic trees were generated from (A) the RdRP coding sequences on the L segment, (B) Glycoprotein (GP) coding sequences on the M segment, and (C) nucleoprotein (NP) sequences on the S segment.Bootstrap values higher than 50 are shown adjacent to the tree branches.The tree is drawn to scale with branch lengths representing the number of substitutions per site.Viral sequences identified from mosquitoes in this study are marked with black circles.The organisms in which the viruses were identified are color-coded with a caption in the box at the right lower corner of the Figures.The viruses considered to be mammalian pathogens or plant pathogens are highlighted in light orange and light green colors, respectively.

Table 1 .
Bunyavirus-positive mosquito species collected in Zambia.*Total RNA sequencing was performed.**Identified only by RNA sequencing analyses.

Table 2 .
Prevalence of the Culex bunyavirus 2 in the Culex quinquefasciatus in Zambia.*The infection rate in the Culex quinuquefasciatus was calculated by the maximum likelihood estimation.**The confidence interval (CI, 95%) was indicated in the range between the upper and lower bounds.