Metatranscriptomic Sequencing Reveals Host Species as an Important Factor Shaping the Mosquito Virome

This study revealed the huge capability of mosquitoes in harboring a rich diversity of RNA viruses, although relevant studies have characterized the intensively unparalleled diversity of RNA viruses previously. Furthermore, our findings showed discernible differences not only in viromic structure between mosquito genera and even between mosquito species within the same genus but also in the genetic diversity and abundance of Wolbachia between different mosquito populations. ABSTRACT Mosquitoes are important vector hosts for numerous viral pathogens and harbor a large number of mosquito-specific viruses as well as human-infecting viruses. Previous studies have mainly focused on the discovery of mosquito viruses, and our understanding of major ecological factors associated with virome structure in mosquitoes remains limited. We utilized metatranscriptomic sequencing to characterize the viromes of five mosquito species sampled across eight locations in Yunnan Province, China. This revealed the presence of 52 viral species, of which 19 were novel, belonging to 15 viral families/clades. Of particular note was Culex hepacivirus 1, clustering within the avian clade of hepaciviruses. Notably, both the viromic diversity and abundance of Aedes genus mosquitoes were significantly higher than those of the Culex genus, while Aedes albopictus mosquitoes harbored a higher diversity than Aedes aegypti mosquitoes. Our findings thus point to discernible differences in viromic structure between mosquito genera and even between mosquito species within the same genus. Importantly, such differences were not attributable to differences in sampling between geographical location. Our study also revealed the ubiquitous presence of the endosymbiont bacterium Wolbachia, with the genetic diversity and abundance also varying between mosquito species. In conclusion, our results suggested that the mosquito host species play an important role in shaping the virome’s structure. IMPORTANCE This study revealed the huge capability of mosquitoes in harboring a rich diversity of RNA viruses, although relevant studies have characterized the intensively unparalleled diversity of RNA viruses previously. Furthermore, our findings showed discernible differences not only in viromic structure between mosquito genera and even between mosquito species within the same genus but also in the genetic diversity and abundance of Wolbachia between different mosquito populations. These findings emphasize the importance of host genetic background in shaping the virome composition of mosquitoes.

which are insect specific (5,6) and are distinct from viruses of medical importance. As such, the presence of disease-causing viruses in mosquitoes is probably the exception rather than the rule. However, how ecological factors impact the composition and diversity of mosquito viromes remains unclear (3,7), limiting our understanding of the ecological drivers of host-switching and spillover events and their public health risks (8). In addition, the pathogenicity of most mosquito-borne viruses is poorly understood, even though some may pose a risk to public health or modulate the transmission of pathogenic viruses (9).
Jinghong City, part of the Xishuangbanna Dai Autonomous Prefecture of Yunnan Province, is located on the southern border of China and geographically adjacent to Laos, Myanmar, and Vietnam, where some tropical viral infectious diseases are endemic (10,11). In recent years, mosquito-associated infectious diseases, such as dengue fever and Japanese encephalitis (12,13), have been imported into Xishuangbanna. In addition, the stable tropical climate in Xishuangbanna facilitates mosquito breeding and production (10), constituting an important driver of infectious disease outbreaks. Indeed, several largescale dengue outbreaks have been reported in Xishuangbanna since the first documented epidemic there in 2013 (12,14,15). Therefore, characterization of the virus spectrum in the key mosquito species (particularly species from the Aedes and Culex genera that are known to harbor pathogenic viruses responsible for epidemics in human and animal populations) and identification of possible associations between mosquito species and virome structures are of vital importance to prevent and control future tropical disease outbreaks here and potentially elsewhere.
Herein, we characterized the total transcriptomes of 56 mosquito pools, comprising 991 mosquitoes from five invertebrate species collected from eight locations in Jinghong, Xishuangbanna. We analyzed the genetic diversity of RNA viruses in these mosquito species and identified the complete coding sequences of 52 RNA viruses, including 19 previously undescribed viruses. We further determined the evolutionary relationships of the novel viruses identified here and revealed an association between mosquito vector species and virome structure by comparing the compositions and structures of the viral communities within different hosts.
For each library, the number of virus species varied from 2 to16, with the exception of one C. quinquefasciatus library from location H, in which no viruses were detected (Fig. 1C). The abundance of each virus varied from 2.00 to 6,825.80 reads mapped per million input reads (RPM) across the pools ( Fig. 1D; Table S2). In comparison, the abundance of the mosquito host, determined by COI gene sequencing, varied between 1. 30 (Table S2). Phylogenetic relationships and viral genome characterization. (i) Double-stranded RNA viruses. We identified 10 double-stranded RNA (dsRNA) viruses belonging to Chrysoviridae (n = 1), Partitiviridae (n = 5), Reoviridae (n = 1), and Totiviridae (n = 3). With the exception of Aedes reo-like virus 1 from Reoviridae, the other nine viruses were related to those previously identified in mosquitoes (Fig. 2).
Two novel ArPLV1 and AePLV2 viruses clustered with uncharacterized partiti-like viruses identified from different mosquito genera ( Fig. 2A), thereby expanding the known host range of partiti-like viruses. The single-gene segments of both viruses contained a single open reading frame (ORF) sharing 83.47% and 61.17% amino acid similarity with their closest relatives, respectively (Table 1). Like other totiviruses, the Culex toti-like virus 1 (CTLV1) identified in the present study possessed an unsegmented genome comprising two major ORFs (Fig. 2B). CTLV1 clustered with unclassified toti-like viruses also from Culex mosquitos and shared 73.20% amino acid similarity with its closest relative (Table 1). Finally, the only novel reovirus identified here, Aedes reo-like virus 1 (AeRLV1), comprised two segments containing three ORFs and was most closely related to Shenzhen reo-like virus 2 from Tyrophagus (Fig. 2D). AeRLV1 shared only 30% amino acid similarity over the conserved RdRp region (Table 1) and formed a distant clade with viruses identified from the class Arachnoidea (Fig. 2D).
(ii) Negative-sense, single-stranded RNA viruses. Twenty-two negative-sense, single-stranded RNA (-ssRNA) viruses with complete coding regions were identified in this study: 7 fell within the order Bunyavirales, 9 fell within the order Mononegavirales, and the remaining 6 belonged to the family Orthomyxoviridae and a related clade (Fig. 3).
Five of the seven Bunyavirales viruses were related to previously described mosquito viruses, with 96 to 100% amino acid identities to their closest relatives, including two Phasmaviridae viruses and three Phenuiviridae viruses (Fig. 3A). Although Culex quinquefasciatus bunyavirus 1 (CQBV1) fell within the Bunyavirales and shared 88.78% amino acid similarity to the most closely related virus (Qingnian mosquito virus [QMV]) (4), both CQBV1 and QMV diverged extremely from previously described viruses and formed a distinct lineage within the Bunyavirales (Fig. 3A), exhibiting less than 23% amino acid similarity over the RdRp protein. Similar to the typical genomic structure of the order Bunyavirales, the two viruses contained three gene segments. Hence, CQBV1 and QMV may represent a putative new family of the Bunyavirales (Fig. 3A). Another novel bunyavirus, Aedes bunya-like virus 1 (ABLV1), identified here also contained three gene segments and fell within the Phenuiviridae (Fig. 3A). ABLV1 clustered with reference strains identified in various mosquitoes, exhibiting 48.54% amino acid similarity to the most closely related virus (Salarivirus) (Fig. 3A).
One novel Armigeres orthomyxo-like virus 1 (ArOLV1) was identified in Armigeres mosquitoes and belonged to an unclassified mosquito-associated clade related to the Orthomyxoviridae (Fig. 3C). ArOLV1 clustered with Usinis virus isolated from Aedes mosquitoes, with an amino acid similarity of 69.27% in the RdRp protein. Although the number of genome segments in the family Orthomyxoviridae ranged from 6 to 8, only 5 genome segments were obtained through sequence similarity from our data, including PB1, PB2, PA, nucleoprotein (N), and the glycoprotein (G) genes (Fig. 3C). (iii) Positive-sense, single-stranded RNA viruses. A total of 20 positive-sense, singlestranded RNA (1ssRNA) viruses were discovered in this study, grouping within the families Tombusviridae (n = 1), Solemoviridae (n = 6), Tymoviridae (n = 2), Flaviviridae (n = 3), Narnaviridae (n = 3), Virgaviridae (n = 2), Permutotetraviridae (n = 1), and the clade Negevlike viruses (n = 2) (Fig. 4). With the exception of Culex tombus-like virus 1 (CTomLV1) within the Tombusviridae, the remaining 1ssRNA viruses identified here were most closely related to mosquito-associated viruses (Fig. 4). CTomLV1 was most closely related to Dansoman virus identified in flies, and the viruses exhibited 48.81% amino acid similarity to each other (Table 1). CTomLV1 had the same genomic structure as Dansoman virus, containing two segments: segment 1 carrying two ORFs of a hypothetical protein and RdRp and segment 2 encoding putative capsid protein (Fig. 4A). We identified two novel viruses of the family Solemoviridae: Aedes sobemo-like virus 1 (ASLV1) and Mosquito sobemo-like virus 1 (MSLV1), which were related to Hubei sobemolike virus 41, previously identified in mosquitoes from China, and Guadeloupe mosquito virus from the Caribbean, respectively (Fig. 4A). ASLV1 and MSLV1 had similar genome structures to their closest relatives and exhibited RdRp amino acid identities of 87.6% and 83.45%, respectively (Table 1). Notably, MSLV1 was identified in both Aedes and Armigeres mosquitoes (Fig. 4A), showing broader host ranges (host-sharing events) across different mosquito genera.
Another host-sharing event was identified in the novel Mosquito narna-like virus 1 (MNLV1), which was also found in both Aedes and Armigeres mosquitoes (Fig. 4B). MNLV1 was closely related to reference strains identified in Aedes, Culex, Coquillettidia, and Ochlerotatus mosquitoes and shared 36.17% amino acid similarity with Ochlerotatus-associated narnalike virus 2. MNLV1 had the same genome structure as its closest relative, with a dual-coding genome structure: two ORFs cover both the sense and antisense genomes, encoding RdRp and a hypothetical protein (Fig. 4B).
We also identified two novel viruses related to the Tymoviridae and Negev-like viruses ( Fig. 4C and D). Culex quinquefasciatus tymo-like virus 1 (CQTLV1) clustered within the clade of mosquito-associated viruses identified in Culex mosquitoes (Fig. 4C), sharing 87.01% amino acid similarity with Culex pseudovishnui tymo-like virus. Similarly, Culex quinquefasciatus negev-like virus 1 (CQNeLV1) and Culex mosquito-associated negev-like viruses formed a clade and showed 90.05% amino acid similarity with Parry's Creek negev-like virus 1 (Fig. 4D). Of particular note was the identification of a novel hepacivirus, termed Culex hepacivirus 1 (CHV1), from a Culex quinquefasciatus mosquito library (Fig. 4G). CHV1 was related to a previously described mosquito-associated virus (Jogalong virus [JgV]) (16), which was identified in Culex annulirostris mosquitoes from Western Australia, sharing 54.62% amino acid sequence similarity with each other. Notably, CHV1 and JgV clustered within the clade of hepaciviruses associated with avian hosts (Fig. 4G), sharing less than 43% amino acid similarity with its closest relative, Bald eagle hepacivirus (Fig. 4G).
Factors affecting the structure and abundance of mosquito viromes. Virus compositions and abundances differed substantially between mosquito species. In general, Aedes mosquitoes (both Ae. aegypti and Ae. albopictus) contained more viruses than Culex mosquitoes (Fig. 1C) and also had higher abundance (Fig. 1D). The highest viral richness was found in Ae. albopictus mosquitoes, with a median of 9.0 (Fig. 5A), followed by Ae. aegypti and C. quinquefasciatus mosquitoes, with medians of 7 and 6, respectively (Fig. 5A). Similarly, both the Shannon (Fig. 5B) and Simpson (Fig. 5C) effective indices were the highest in Ae. albopictus, followed by Ae. aegypti, and the lowest values of the two indices were from C. quinquefasciatus (Fig. 5).
Of the 52 viral species discovered, 5 were shared between Aedes and Culex mosquitoes ( Fig. 5D; Table S2), and only 2 were shared between Aedes and Armigeres species (Fig. 5D). No viruses were shared between Culex and Armigeres mosquitoes (Fig. 5D). Notably, Ae. aegypti shared 14 viruses with Ae. albopictus, while C. quinquefasciatus shared one virus with other Culex species (Lutzia halifaxii), although only two viruses were discovered in L. halifaxii (Fig. 5D). The clustering of the libraries with similar viromic composition and abundance was further described using b-diversity analysis. Principal-coordinate analysis (PCoA) plots based on the Bray-Curtis distance matrices revealed that samples from the same mosquito species clustered together (adonis; R 2 = 0.67, P = 0.0002) (Fig. 6A), demonstrating that mosquito species affect virome structure.
The possible association between sampling locations and viral composition and abundance was examined further (Fig. S1) by measuring observed richness (Fig. S1A), the Shannon index (Fig. S1B), and the Simpson index (Fig. S1C). Three locations (D, E, and F) contained only one or two libraries with limited virome diversity. Only location H showed a significant difference from locations B and C, suggesting a potential association between the geographic location and the viromic structure of mosquitoes. In addition, PCoA plots based on the Bray-Curtis distance matrices revealed that samples from different locations did not form distinct clusters (Fig. S1D). Hence, there is no strong evidence for geographic structure in mosquito viromes in this study, and the mosquito collection sites were close to each other in this study.
We further constructed a co-occurrence network among mosquito species based on significant positive correlations (Spearman's r . 0.6; P , 0.05) (Fig. 6B). The network was derived from the relative abundance of each virus, comprising 56 nodes (56 mosquito libraries) and 372 edges. Based on the modularity class, the entire network could be parsed into three major modules, corresponding to the three major mosquito species (Fig. 6B). Notably, nodes were inclined to interact more with nodes within the same module than with nodes of other modules. The co-occurrence patterns clearly illustrated the correlations between the mosquito species and viral abundance; however, the network among sampling locations showed no obvious correlations (Fig. S2).
Wolbachia diversity and abundance. As an endosymbiont bacterium, Wolbachia was detected in all libraries in this study, with relatively high abundance (11.72 to 2,899.50 RPM) (Fig. 1D; Table S2). The alignment of Wolbachia 16S rRNA sequences showed 87.99 to 100% nucleotide identities to each other, and phylogenetic analyses revealed two major lineages (Fig. 7A): the first comprised all Wolbachia sequences obtained from Ae. aegypti, while the other contained Wolbachia sequences from all five mosquito species. The phylogenetic tree indicated that Wolbachia sequences from Ae. aegypti had higher genetic diversity, while Wolbachia sequences from Ae. albopictus and C. quinquefasciatus were more genetically homogeneous. Further comparisons of the abundance (RPM) of Wolbachia gene sequences across libraries revealed significant differences between Ae. aegypti, Ae. albopictus, and C. quinquefasciatus (Fig. 7B). The highest abundance of Wolbachia sequences was observed   in C. quinquefasciatus, with a median of 1,599.5 RPM, followed by Ae. albopictus mosquitoes (median of 977.8), and the lowest abundance was in Ae. aegypti mosquitoes, with a median of 25.5 (Fig. 7B).

DISCUSSION
We present a comprehensive description of the viromes of 991 mosquitoes collected from eight locations in Yunnan Province in southwestern China. Although previous metagenomic studies have revealed numerous novel and highly divergent RNA viruses in mosquitoes, analysis of the transcriptomes of the five mosquito species in the present study has led to the discovery of 52 RNA viruses belonging to 15 viral families and unclassified clades. Notably, 19 viruses were novel, with half sharing 30 to 70% amino acid similarity to their most closely related viruses. One of the most notable discoveries was Culex quinquefasciatus bunyavirus 1 (CQBV1), which represented a putative new virus family within the order Bunyavirales together with previously described QMV (4). These results underscore the capacity of mosquitoes to harbor a wide diversity of RNA viruses, highlighting the necessity for constant surveillance of potential viral pathogens in these arthropod vectors.
Another key result from our study was that within the Aedes mosquito vector species, both Ae. aegypti and Ae. albopictus harbored significantly greater virus diversity than Culex mosquitoes (Culex quinquefasciatus). This is in contrast to the observation in a previous report that Culex species harbored more viruses at high abundance than Aedes mosquitoes (17). The underlying explanation for these contrary results could be due to uneven sample sizes of each mosquito genus. Our results also revealed pronounced differences between the virome structures of Ae. aegypti and Ae. albopictus mosquitoes: even though there is considerable overlap in the viruses carried by these two Aedes mosquitoes, the latter had both higher diversity and viral abundance. These findings were also supported by further statistical and multidimensional scaling analyses and are consistent with prior evidence that different Aedes mosquitoes can have significantly different virome compositions (17,18). Hence, the ecology of mosquito viruses was driven, at least in part, by host taxon, consistent with the predicted narrow host range of insect-specific viruses (6).
In contrast, according to the aand b-diversity analyses, the virome structures were relatively homogeneous across different locations. Indeed, some mosquito-specific viruses exhibiting high similarity have been described from different continents (19), indicating that viruses can be transmitted to wide geographical areas through mosquito populations. As all of the sampling sites in the present study were geographically close, larger-scale samplings covering different ecological niches are clearly required for further investigations with respect to the correlation between geographic location and virome structure.
The viruses discovered here expand the host ranges of several mosquito viruses to include additional mosquito species and even genera. For example, Shuangao chryso-like virus 1 YN2018 (ShCLV1), Culex phasma-like virus YN2018 (CPLV), and Guangzhou sobemolike virus YN2018 (GSLV) were identified in both Aedes and Culex mosquitoes, while Mosquito narna-like virus 1 (MNLV1) was present in both Aedes and Armigeres mosquitoes with high abundance, suggesting host sharing and the intergeneric transmission of these viruses. However, those viruses that clustered with viruses associated with fungi rather than mosquitoes or arthropods might have been derived from other eukaryotic organisms present in the mosquito microbiome or from fungal infections of the mosquito cuticle.
According to the WHO (20), the most prevalent viral infections are primarily transmitted by Ae. aegypti, Ae. albopictus, and C. quinquefasciatus, including Zika virus fever, dengue, yellow fever, Japanese encephalitis, and West Nile fever. However, surprisingly, none of these known viral pathogens were identified in this study. Interestingly, one virus, CHV, as well as the previously documented JgV (16), clustered with vertebrate-associated hepaciviruses. Phylogenetic analyses of the Flaviviridae suggested that both CHV and JgV were likely associated with avian hosts, rather than the mosquito itself. Indeed, through targeting two of the vertebrate mitochondrial genes (COI and Cytb), JgV was suspected as a contamination from a blood meal taken from a bird host (16). However, we did not find vertebrate mitochondrial genes from the sequencing data. Given the high divergence of these mosquito hepaciviruses, it will be important to investigate their true natural host, particularly as this may provide valuable information on their evolutionary history. Although analysis of host genes from mosquito sequencing data is highly suggestive of the natural host, viruses detected in the blood of vertebrate species would provide convincing evidence of the real source of the viruses. To this end, more expansive surveillance is required, including a larger collection of mosquitoes and blood samples from vertebrate animals present in the same location.
Wolbachia has been documented to provide resistance to the infection with some viruses, such as dengue virus and Zika virus in mosquitoes (21,22), and hence has been suggested as a potential tool for vector-borne disease control. However, not all Wolbachia strains have clear effects in inhibiting virus replication, and Wolbachia infection does not protect against all viruses (23,24). Despite this, little is known about how the ecology of hosts impacts Wolbachia diversity. We found Wolbachia in all libraries, once again indicating the prevalence of Wolbachia in mosquitoes in China (25). However, there was a large discrepancy in the genetic diversity and abundance of Wolbachia between different mosquito populations. Specifically, Wolbachia in Ae. aegypti had high diversity but low abundance, while the converse was seen (low diversity/high abundance) in both Ae. albopictus and C. quinquefasciatus. These results suggested that mosquito species might also play an important role in Wolbachia composition.
We do not believe that the Wolbachia sequences identified in Ae. Aegypti mosquitoes result from contamination as the Ae. aegypti libraries were sequenced on different lanes and sequences within the same lane did not share 100% nucleotide identity. In addition, Wolbachia-infected Ae. Albopictus mosquitoes, not Ae. Aegypti, are being released in a few small independent islands in China (26), and our samples were not collected from these sites. However, Ae. aegypti was not thought to naturally harbor Wolbachia (27) until recently in some Southeast Asian countries (such as Malaysia, India, the Philippines, and Thailand) as well as the United States (28,29). None of these studies provide robust evidence that Ae. aegypti harbors natural Wolbachia infections. The presence of natural Wolbachia infections may interfere with compatibility patterns between Ae. aegypti mosquitoes and some Wolbachia strains (27). Confirmation of a natural infection in these mosquitoes will require significant additional experimental work.
In conclusion, we studied the viral diversity of five mosquito species sampled in different locations in Yunnan Province and highlighted the capacity of mosquitoes to harbor a rich diversity of RNA viruses. The viral compositions varied mainly between different mosquitoes, suggesting host species represents an important factor shaping the virome composition of mosquitoes.

MATERIALS AND METHODS
Sample collection. From October to December 2018, a total of 991 adult mosquitoes were collected using light-traps from eight locations in Jinghong City, Xishuangbanna, Yunnan Province, China (Fig. 1). Mosquito species were initially identified morphologically by experienced field biologists and further confirmed based on sequences of the cytochrome c oxidase subunit 1 mitochondrial gene (COI). All samples were divided into 56 pools by mosquito species and geographic location and were transported to the laboratory on dry ice.
Metatranscriptomic sequencing. Mosquitoes were rinsed three times using RNA-and DNA-free phosphate-buffered saline (PBS) solution (Gibco) before homogenization with steel beads in PBS solution. Total RNA was extracted using TRIzol reagent (Invitrogen). RNA quantity and quality were checked using the Agilent 2100 Bioanalyzer (Agilent). Host rRNA was removed using the MGIEasy rRNA depletion kit according to the manufacturer's instructions, and sequencing libraries were constructed using the MGIEasy mRNA library prep kit. Pairedend (100-bp) sequencing of each RNA library was performed on the BGISEQ-500RS sequencing platform (BGI).
Data analysis and virus discovery. A quality assessment of the raw sequencing reads was conducted using the Fastp v.0.19. (30) and Trimmomatic (31) programs, before de novo assembly using the Trinity program (32). The assembled contigs were then compared against the nonredundant nucleotide (nt) and protein (nr) databases downloaded from NCBI using blastn (33) and Diamond blastx (34), with cut-off E values of 1 Â 10 210 and 1 Â 10 25 , respectively. All potential viral contigs were identified and merged into longer viral contigs using Geneious Prime (35). False-positive results due to cross-contamination and index hopping during sequencing were excluded as previously described (36). The relative abundance of the viruses identified was determined by mapping the reads back to the assembled contigs using Bowtie2 v.2.3.3.1 (37) and calculated as the number of reads mapped per million input reads (RPM) using the formula "total mapped reads/total reads Â 1 million." Bowtie2 was used to align the reads to each novel virus genome, and SAMtools (38) was used to compute the percentage of reads mapped and coverage depth. Novel viruses were defined employing the previously Role of the Host in Shaping the Mosquito Virome Microbiology Spectrum defined criterion such that the translated protein sequence shared less than 90% amino acid similarity in the RNA-dependent RNA polymerase (RdRp) to any previously described viruses (39). Phylogenetic analyses. RdRp sequences of the viruses identified from this study were then aligned with their corresponding homologs in reference viruses using the MAFFT v.7.407 program (40) employing the E-INS-I algorithm, followed by the removal of ambiguously aligned regions using TrimAl v.1.4 (41). Phylogenetic trees were constructed using the maximum likelihood method implemented in IQ-TREE v.1.6.12 (42), employing the best-fit substitution models identified by IQ-TREE.
Ecological dynamics analysis and network analysis. Statistical analyses of viral genetic diversity and abundance were performed using the t test or Wilcoxon test based on the results of a normal distribution test (Shapiro-Wilk test) in the ggpubr package and were plotted using the ggplot2 package in RStudio v.4.1.2. The observed richness, Shannon index, and Simpson index (i.e., a diversity) were estimated for each library using modified Rhea script sets (43) and compared between different mosquitoes using the Kruskal-Wallis rank sum test. Principal-coordinate analysis (PCoA) (i.e., b diversity) was performed based on the Bray-Curtis dissimilarity matrix using the Vegan package (44). A correlation between two items was considered statistically robust if the Spearman's correlation coefficient (r ) was .0.6 and the P value was ,0.05, with the P value adjusted with a multiple-testing correction using the Benjamini-Hochberg method (45). All pairwise Spearman's rank correlations between the viral members were calculated using the psych package in RStudio v.4.1.2. Network visualization was conducted on the interactive platform of Gephi (46).
Identification of Wolbachia bacteria. Bacteria were initially identified in the metatranscriptomic data using MetaPhlAn2 (47). The 16S rRNA gene was used to conduct phylogenetic analyses and similarity comparisons for Wolbachia. To estimate their relative abundance, sequence reads were mapped to the complete reference genomes (CP031221 [Wolbachia pipientis wAlbB chromosome]) from which the RPM was calculated.
Ethics statement. This study was performed in accordance with the institutional and national guidelines for the care and handling of the animals.
Data availability. All sequence reads generated in this study have been uploaded into the NCBI Sequence Read Achieve (SRA) database under BioProject accession no. PRJNA911492. All novel and known virus genome sequences generated in this study have been deposited in NCBI/GenBank under accession no. OQ067620 to OQ067711.

SUPPLEMENTAL MATERIAL
Supplemental material is available online only. SUPPLEMENTAL FILE 1, PDF file, 1.1 MB. SUPPLEMENTAL FILE 2, XLSX file, 0.04 MB. We declare no conflict of interest.