Meta-Transcriptomic Comparison of the RNA Viromes of the Mosquito Vectors Culex pipiens and Culex torrentium in Northern Europe

Mosquitoes harbor an extensive diversity of ‘insect-specific’ RNA viruses in addition to those important to human and animal health. However, because most studies of the mosquito virome have been conducted at lower latitudes, little is known about the diversity and evolutionary history of RNA viruses sampled from mosquitoes in northerly regions. Here, we compared the RNA virome of two common northern mosquito species, Culex pipiens and Culex torrentium, collected in south-central Sweden. Following bulk RNA-sequencing (meta-transcriptomics) of 12 libraries, comprising 120 specimens of Cx. pipiens and 150 specimens of Cx. torrentium, we identified 40 viruses (representing 14 virus families) of which 28 were novel based on phylogenetic analysis of the RNA-dependent RNA polymerase (RdRp) protein. Hence, we documented similar levels of virome diversity as in mosquitoes sampled from the more biodiverse lower latitudes. Many viruses were also related to those sampled on other continents, indicative of a widespread global movement and/or long host–virus co-evolution. Although the two mosquito species investigated have overlapping geographical distributions and share many viruses, several viruses were only found at a specific location at this scale of sampling, such that local habitat and geography may play an important role in shaping viral diversity in Culex mosquitoes.


Introduction
The mosquito genus Culex (Diptera; Culicidae) comprises more than a thousand species, with representatives found globally [1]. Culex species are vectors of a number of important pathogens including West Nile virus (WNV) (Flaviviridae), Japanese encephalitis virus (JEV) (Flaviviridae), and Sindbis virus (SINV) (Togaviridae), as well as a variety of nematodes [1][2][3]. One of the most widespread Culex species is the Northern House mosquito, Cx. pipiens, which is distributed across the northern hemisphere. In Europe and the Middle East, it occurs together with Cx. torrentium, another Culex species with females and larvae that are morphologically identical to Cx. pipiens. These two species have overlapping distributions and share larval habitats. However, Cx. torrentium dominates in Northern Europe while Cx. pipiens is more abundant in the south [4]. Both species are vectors for a number of bird-associated viruses that can cause disease in Europe, such as WNV, that may cause a febrile disease with encephalitis, and SINV that may result in long-lasting arthritis [2,5]. Cx. pipiens is one of the most common WNV vectors in both Southern Europe and North America, while Cx. torrentium is the main vector of SINV in Northern Europe [2,6]. Infections with these pathogenic viruses occur in late summer when viral prevalence increases in passerine birds, the vertebrate hosts of both of these viruses [7,8]. Despite their importance as vectors, little is known about the detailed biology of Cx. pipiens and Cx. torrentium due to the difficulties in species identification, which can only be reliably achieved through molecular means. Much of the biology of these species, such as their larval habitat and feeding preferences, is considered similar. However, one important difference between the two species is that while Cx. pipiens harbors a high prevalence of the intracellular bacteria Wolbachia pipientis, it is seemingly absent in Cx. torrentium [9].
In recent years, studies utilizing RNA-sequencing (RNA-Seq, or 'meta-transcriptomics') have revealed an enormous RNA virus diversity in both vertebrates and invertebrates [10,11]. Mosquitoes are of particular interest as many are well-known vectors of pathogenic viruses. Importantly, these pathogenic viruses represent only a fraction of the total virome in the mosquito species investigated. Indeed, mosquitoes clearly carry a large number of newly described and divergent arthropod-specific viruses, with representatives from many genetically diverse virus families and orders, such as the Flaviviridae, Togaviridae, and the Bunyavirales [12][13][14][15][16]. However, most studies have been conducted on latitudes below 55 • , such that there is a marked lack of data of the mosquito viral diversity present in northern temperate regions where the composition of mosquito species as well as environmental parameters differ significantly from lower latitudes, and where human populations are at high density. In addition, for many life forms, biodiversity increases towards the equator [17], and the species richness of mosquitoes is greater in tropical regions than temperate regions [18]. A central aim of the current study was therefore to investigate whether viral diversity co-varies in the same manner. Given that Cx. pipiens and Cx. torrentium are two common Culex species in Northern and Central Europe, and known vectors of SINV and WNV, they were chosen for RNA virome investigation and comparison by RNA-Seq.

Mosquito Collection
Mosquitoes were collected from two regions in Sweden: (i) from floodplains of the Dalälven River in central Sweden (60.2888; 16.8938) in 2006, 2009, and 2011; and (ii) around the city of Kristianstad, in southern Sweden (56.0387; 14.1438) in 2006 and 2007. Collections were performed using Centers for Disease Control and Prevention-light traps baited with carbon dioxide, and catches were sorted and identified to species on a chilled table, using keys by Becker et al. [19]. In total, legs from 270 Cx. pipiens/torrentium mosquitoes were removed to enable molecular identification to species [20]. Bodies were homogenized in phosphate-buffered saline buffer supplemented with 20% fetal calf serum and antibiotics and stored at -80 • C until further processing.

Sample Processing and Sequencing
Total RNA was extracted from 12 pools from the homogenate of individual Cx. torrentium (n = 150) and Cx. pipiens mosquitoes (n = 120) (Table S1), using the RNeasy ® Plus Universal kit (Qiagen, Hilden, Germany) following the manufacturer's instructions. Three pools, L1 and L2 for Cx. torrentium and L3 for Cx. pipiens, respectively, were previously shown to be positive for SINV and were included as a reference for RNA viral diversity in the presence of a pathogen. The extracted RNA was subsequently DNased and purified using the NucleoSpin RNA Clean-up XS kit (Macherey-Nagel, Düren, Germany). Prior to library construction, ribosomal RNA (rRNA) was depleted from the purified total RNA using the Ribo-Zero Gold (human-mouse-rat) kit (Illumina, San Diego, CA, USA) following the manufacturer's instructions. Sequencing libraries were then constructed for all rRNA-depleted RNA-samples using the TruSeq total RNA library reparation protocol (Illumina). All libraries were sequenced on a single lane (paired-end, 150 bp read-length) on an Illumina HiSeq X10 platform. Library preparation and sequencing was carried out by the Beijing Genomics Institute, Hong Kong. All 12 libraries were quality trimmed with Trimmomatic v.0.36 [21], using default settings for paired-end sequence data, and then assembled de novo using Trinity v.2.5.4 [22], employing the default settings with read normalization.

Identification of Viruses and Wolbachia Bacteria
Trinity assemblies were screened against the complete non-redundant NCBI GenBank nucleotide (nt) and protein (nr) databases using blastn, primarily to identify closely related RNA viruses and false-positive host-derived hits, as well as a diamond [23] blastx analysis primarily to identify divergent RNA viruses, with cut-off e-values of 1 × 10 −5 in both cases. To determine whether some of the assemblies represent potential endogenous viral elements (EVEs), all virus viral assemblies were blasted against the Culex quinquefasciatus reference genome (GCA_000209185.1). Assemblies identified as RNA viruses were screened against the NCBI Conserved Doman Database with an expected value threshold of 1 × 10 −3 to identify viral sequence motifs. The mitochondrial cytochrome c oxidase I (COX1) gene, mined from the sequence data, and all contigs with RdRp-motifs was mapped back to all quality trimmed libraries to estimate abundance using Bowtie2 [24], employing the default local setting. A virus was considered to be in high abundance if: (i) it represented >0.1% of total ribosomal-depleted RNA reads in the library, and (ii) if the abundance was higher to that of the abundant host COX1 gene [12,25]. Such high abundance viruses were tentatively assumed to be mosquito associated. Hits below the level of cross-library contamination due to index-hopping, measured as 0.1% of the most abundant library for the respective virus species or less than 1 read per million mapped to a specific virus contig, was considered negative (colored grey in Tables 1 and 2, respectively). To investigate the presence of Wolbachia bacteria in the libraries, published sequences of the Wolbachia Cx. pipiens wsp surface protein gene (DQ900650.1) and the mitochondrial COX1 gene (AM999887.1) were mapped backed against all libraries using the above criteria for abundance and presence/absence.

Inference of Virus Evolutionary History and Host Associations
The evolutionary (i.e., phylogenetic) history of the viruses discovered were inferred by aligning protein translated open reading frames with representative sequences from the Alphaviridae, (Order) Bunyavirales, Endornaviridae, Luteoviridae, (Order) Mononegavirales, Nido-like viruses, (Order) Orthomyxovirales, Partitiviridae, Picornaviridae, Qin-like viruses, Reoviridae, Totiviridae, Tymoviridae and Virgaviridae and Negev-like viruses. All RdRp amino acid sequence alignments were performed using the E-INS-i algorithm in Mafft [26]. Poorly aligned regions, in which amino acid positional homology could not be confirmed, were then removed from the alignments using TrimAl utilizing the 'strict' settings. Finally, phylogenetic trees were computed with a maximum likelihood approach as implemented in PhyML [27] employing the LG+Γ model of amino acid substitution, Sub-tree Pruning and Re-grafting branch-swapping and the approximate likelihood ratio test (aLRT) with the Shimodaira-Hasegawa-like procedure used to assess branch support. The resultant phylogenetic trees were edited and visualized with FigTree v.1.4.2.
To help assess whether the novel viruses discovered were likely mosquito associated, that is, to distinguish those that actively replicate in the host from those present in diet mosquito or a co-infecting micro-organism, we considered four factors: (i) the abundance of viral contigs per total number of reads in a library (i.e., >0.01% was considered abundant); (ii) if the abundance of the virus was higher in relation to the host COX1 gene; (iii) presence in the individual sequencing libraries (i.e., present = yes, not present = no); and (iv) clear phylogenetic clustering with other mosquito-derived viruses (i.e., clustered/did not cluster with other mosquito associated viruses = yes/no). A mosquito association was tentatively assigned if a particular virus met two or more of the four criteria. The raw sequence data generated here have been deposited in the NCBI Sequence Read Archive (BioProject: PRJNA516782) and all viral contigs have been deposited on NCBI GenBank (accession numbers: MK440619-MK440659).

RNA Virome Characterization
We characterized the RNA viral transcriptome of two mosquito species, Cx. pipiens and Cx. torrentium, collected from central and southern Sweden (Table S1). After high-throughput sequencing, a total of 569,518,520 (range 34,150,856-62,936,342) 150bp reads were produced from 12 ribosomal RNA-depleted sequence libraries that were assembled into 153,583 (4333-33,893) contigs. From all the contigs assembled, we identified 40 that contained RdRp sequence motifs and hence indicative of viruses, belonging to 14 different viral families/orders: Alphaviridae, Bunyavirales, Endornaviridae, Luteoviridae, Mononegavirales, Nidovirales, Orthomyxoviridae, Partitiviridae, Picornaviridae, Reoviridae, Totiviridae, Tymoviridae, and representatives from the divergent Virgaviridae, Negeviridae, and Qin viruses. Following a similarity search of all virus sequences against a Cx. quinquefasciatus reference genome we found no evidence that any of the discovered viruses were derived from the mosquito host genome (i.e., present as EVEs). For each viral family/order, between one and five virus species were identified and in total 28 novel RNA virus species were discovered here, which were named based on geographical location.

Virome Comparison between Mosquito Species and Geographical Regions
Both the composition and abundance of the virus species and families observed seemingly differed between the two mosquito species (Figure 2, Table 2). Of the 40 newly discovered virus species, most were found in Cx. pipiens which harbored 34 species: 23 of these are newly described in Cx. pipiens and 11 have been described previously. Sixteen of these 34 virus species were unique to Cx. pipiens and hence not present in Cx. torrentium. Similarly, 24 of the 40 virus species were discovered in Cx. torrentium: 18 of these were newly described in Cx. torrentium and six have been described previously. Six viruses found in Cx. torrentium were not present in Cx. pipiens ( Figure 3, Table 2).
We next analyzed potential host relationships by comparing the abundance (total abundance, of which >0.01% was considered abundant, and in relation to the host COX1 gene), presence across multiple libraries, and phylogenetic relationship to other viruses (

Virome Comparison between Mosquito Species and Geographical Regions
Both the composition and abundance of the virus species and families observed seemingly differed between the two mosquito species (Figure 2, Table 2). Of the 40 newly discovered virus species, most were found in Cx. pipiens which harbored 34 species: 23 of these are newly described in Cx. pipiens and 11 have been described previously. Sixteen of these 34 virus species were unique to Cx. pipiens and hence not present in Cx. torrentium. Similarly, 24 of the 40 virus species were discovered in Cx. torrentium: 18 of these were newly described in Cx. torrentium and six have been described previously. Six viruses found in Cx. torrentium were not present in Cx. pipiens ( Figure 3, Table 2).
We next analyzed potential host relationships by comparing the abundance (total abundance, of which >0.01% was considered abundant, and in relation to the host COX1 gene), presence across multiple libraries, and phylogenetic relationship to other viruses (Table 3). If a particular virus met two or more of the four criteria, it was tentatively considered as mosquito associated. These data suggest that 16 of the 40 viruses were likely mosquito associated, of which one and two were unique to Cx. pipiens and Cx. torrentium, respectively (Figures 4-6). The host association was unclear in the remaining viruses (for example, they could be associated with micro-organisms co-infecting the mosquitoes) and hence could not be safely assumed to infect mosquitoes. For example, Ahus virus (Totiviridae) was at low abundance, was not present in several libraries, and clustered with viruses  (Table 3), its closest relative ( Figure 5) was a soybean leaf-associated virus [28] which means that its proposed mosquito association is uncertain and clearly needs to be investigated further. Conversely, Culex mononega-like virus 2 ( Figure 5) was abundant, present in several libraries and clustered with other mosquito viruses, suggesting that it is mosquito associated. All potential mosquito host association data is summarized in Table 3.
Viruses 2019, 11, x FOR PEER REVIEW 2 of 22 two or more of the four criteria, it was tentatively considered as mosquito associated. These data suggest that 16 of the 40 viruses were likely mosquito associated, of which one and two were unique to Cx. pipiens and Cx. torrentium, respectively (Figure 4-6). The host association was unclear in the remaining viruses (for example, they could be associated with micro-organisms co-infecting the mosquitoes) and hence could not be safely assumed to infect mosquitoes. For example, Ahus virus (Totiviridae) was at low abundance, was not present in several libraries, and clustered with viruses derived from various environmental samples, suggesting that it is not mosquito associated. Similarly, although Gysinge virus (Mononegavirales) was abundant and present in several libraries (Table 3), its closest relative ( Figure 5) was a soybean leaf-associated virus [28] which means that its proposed mosquito association is uncertain and clearly needs to be investigated further. Conversely, Culex mononega-like virus 2 ( Figure 5) was abundant, present in several libraries and clustered with other mosquito viruses, suggesting that it is mosquito associated. All potential mosquito host association data is summarized in Table 3.
Cx. torrentium Cx. pipiens   Cx. pipiens -Kristianstad # unique viruses = 5 Notably, Cx. torrentium carried four viruses of markedly higher abundance compared to Cx. pipiens: (i) Nam Dinh virus (303,145 RPM, or 42% of all viral reads and more than 30% of all (non rRNA) reads, respectively, in library L2); (ii) Biggie virus (35,063 RPM, or 4.6% of all viral reads and 3.5% of all reads, respectively, in library L12); as well as two newly identified viruses, (iii) Valmbacken virus (52,264 RPM, or 27% of all viral reads and 3.5% of all reads, respectively, in library L12); and (iv) Jotan virus (41,738 RPM, or 56% of all viral reads and 4.2% of all reads, respectively, in library L11) ( Table 2, Table  S2). Cx. pipiens was characterized by a slightly higher abundance of orthomyxo-like and luteoviruses compared to Cx. torrentium (Figure 2), although in both mosquito species the most abundant virus was the Nam Dinh virus that reached 10,006 RPM (or 73% of all viral reads and 1% of all reads, respectively) in library L5. Table 3. Indication of host associations for the viruses discovered here. Likely host association was assessed using (i) the abundance of viral contig per total amount of reads in a library, (ii) virus abundance in relation to the host COX1 gene, (iii) presence across libraries, and (iv) phylogenetic clustering with other mosquito-derived viruses. If a particular virus met two or more of the four criteria, it was considered as a mosquito-associated virus. We next compared the virome composition in the sampled mosquitoes between Kristianstad in the south and the floodplains of the Dalälven River situated roughly 600 km further north. In the case of Cx. pipiens this analysis revealed a total of 20 virus species from Kristianstad, 12 of which were unique to Cx. pipiens and five detected in Kristianstad only, all of which were unique to Cx. pipiens: Asum virus (Bunyaviridae), Eskilstorp virus (Chrysoviridae), Kristianstad virus (Bunyaviridae), Rinkaby virus (Virga-Negev virus), and Vittskovle virus (Qinvirus). A total of 28 viruses were found in Cx. pipiens from Dalälven. Eleven of these were unique to Cx. pipiens and four were unique to Cx. pipiens from Dalälven: Salari virus (Bunyavirales), Sonnbo virus (Partitiviridae), Culex mononega-like virus 1 (Mononegavirales), and Berrek virus (Luteoviridae) ( Table 2, Figures 2 and 4-6). A similar relationship was found for Cx. torrentium. In the case of Dalälven, 24 viruses were found in Cx. torrentium, of which 18 were shared with Cx. pipiens and six of which were unique to Cx. torrentium ( Table 2, Figures 2 and 4-6). Hence, the majority of the mosquito viruses identified here were shared both between species and geographical regions, even though only 30 specimens of Cx. pipiens were available from Kristianstad. In contrast, we found several virus species that were unique to a specific location, which could be indicative of virome differentiation at a local geographic scale, although this will need to be confirmed with additional sampling.

Evolutionary History and Host Associations of the Discovered RNA Viruses
Our phylogenetic analyses of the viruses newly identified here showed that several were closely related to previously identified viruses, and that many form clusters with mosquito associated and/or Culex associated viruses within particular viral families, such as the Merida virus and Gysinge virus (Mononegavirales), and Tarnsjo virus (Tymovirales) (Figures 4-6). In contrast, other novel viruses clustered with those neither associated with mosquitoes nor other arthropods: that they are distinguished by long branches suggests that they might infect diverse host taxa.

Positive-Sense RNA Viruses
We identified 16 positive-sense RNA viruses, of which 12 were likely novel. The majority fell within the Hepe-Virga-Endorna-Tymo-like virus complex (n = 8), whereas the others fell within Nidovirales (n = 1), Luteoviridae (n = 4), Picornavirales (n = 2), and Togaviridae (n = 1), respectively (Figure 4). The viruses discovered contain those that are closely related to other mosquito-associated viruses, such as the highly abundant Nam Dinh virus (Nidovirales), as well as those without clear host associations. For example, Biggie virus clusters in a distinct group of Biggie viruses (Virga/Endorna-viridae) sampled from other Culex mosquitoes [15]. We also identified several novel and divergent viruses in the Endornaviridae-specifically the Kerstinbo virus and Hallsjon virus-that do not cluster with other arthropod-associated viruses (Figure 4). Similarly, within the Tymoviridae we detected two variants of a Culex-associated virus, Tarnsjo virus, that are closely related to a Culex associated Tymoviridae-like virus [15].
We identified four viruses within the Luteoviridae: Culex associated luteo-like virus, as well as the novel Berrek, Fagle, and Marma viruses. Culex associated luteo-like virus was previously found in a pool of Culex sp. mosquitoes from North America [15]. Both the newly discovered Berrek virus and Marma virus grouped with other luteoviruses found in mosquitoes (Figure 4), but only the Marma virus was abundant, suggesting that it is Culex associated ( Table 2, Table 3).
Two novel picornaviruses were also identified. The abundant Rinkaby virus clustered with Yongsan iflavirus 1 virus, sampled from Cx. pipiens mosquitoes from South Korea and was therefore considered a bona fide Culex associated picornavirus. Although the Ista virus did not cluster with viruses derived from mosquitoes, its high abundance and presence in all libraries (Figure 4, Table 3) suggest that it is also Culex associated.
Finally, four of our libraries-L1 and L2 for Cx. torrentium and L3 and L6 for Cx. pipiens-contained reads for SINV. Importantly, whereas the presence of SINV could be confirmed with PCR in library L1, L2 and L3, it was not PCR confirmed in L6 such that contamination cannot be excluded in this case.

Negative-Sense RNA Viruses
In total, we identified 16 negative-sense RNA viruses, including nine novel viruses: Bunyavirales (n = 5), Mononegavirales (n = 5), Qin-like viruses (n = 4), and Orthomyxoviridae (n = 2) ( Figure 5). As was the case for the positive-sense RNA viruses, some of these viruses have been identified previously and clustered with viruses found in mosquitoes of the same genera, including Salari virus (Bunyavirales) and a number of novel viruses such as Anjon virus (Bunyavirales).
Within the order Mononegavirales, Merida virus, Culex mononega-like virus 2, Culex mosquito virus 4, and Culex mononega-like virus 1 have previously been described in mosquitoes [12,15,29]. Although abundant, the novel Gysinge virus did not cluster with any mosquito sequences ( Figure 5), so its true host is uncertain.
The Qinviruses are a newly described and highly divergent group of RNA viruses [10]. We identified four novel Qin-like viruses: Nackenback virus, Gran virus, Vinslov virus, and Vittskovle virus. The latter three are more closely related to the Hubei qinvirus-like virus 2 previously found in a pool of different arthropod species [10]. Nackenback virus was found to share a more recent common ancestor with the Wilkie Qin-like virus previously found in Aedes and Culex mosquitoes in Australia [12]. Although Qin-like was most closely related to fungal viruses [12], it is notable that Nackenback virus was found in both Cx. pipiens and Cx. torrentium libraries and was also more abundant than host non-RNA in the Cx. torrentium libraries ( Table 2, Table 3). Hence, this virus may be truly mosquito associated. We also detected two orthomyxoviruses, Wuhan Mosquito Virus 6 and Wuhan Mosquito Virus 4, both of which have previously been found in pools of Culex mosquitoes and are known to be mosquito associated [12].

Double-Stranded RNA Viruses
We identified eight double-stranded RNA viruses in our Swedish mosquitoes, all of which were novel: Partitiviridae (n = 2), Reoviridae (n = 1), and Toti/Chrysoviridae (n = 5). For the family Reoviridae, Valmbacken virus clustered with Aedes camptorhynchus reo-like virus, previously discovered in mosquitoes [12]. Valmbacken virus was also abundant and found in all libraries and is therefore most likely a Culex associated reovirus ( Table 2, Table 3).
In comparison, four of the five novel totiviruses clustered with other mosquito-associated totiviruses ( Figure 6), but only two (Lindangsbacken virus and Eskilstorp virus) were also abundant. The fifth totivirus, Ahus virus, was highly divergent, had low abundance, and clustered distantly with potentially protist originating viruses ( Figure 6). Thus, the host association of Ahus virus remains uncertain.
The two novel partiti-like viruses, Vivastbo virus and Sonnbo virus, did not cluster with any viruses sequenced from mosquitoes, but rather grouped with viruses originating from various arthropod hosts ( Figure 6). However, the relatively high abundance levels of the Vivastbo virus suggest that it may be associated with mosquitoes ( Table 2, Table 3).
Finally, all sequencing libraries generated here were negative for Wolbachia as assessed by mapping against the COX1 and Wolbachia surface protein (WSP) genes of Wolbachia pipientis. Although it is not possible to completely exclude the presence of other Wolbachia variants, our results suggest that differential presence/absence of Wolbachia has not affected the observed patterns of viral diversity and abundance.

Discussion
Through total RNA-sequencing of 270 Culex mosquitoes collected in Sweden we identified 40 viruses, including 28 that are novel. A virome comparison between the two vector species Cx. pipiens and Cx. torrentium revealed that although these mosquitoes are from the same genus and have overlapping geographical distribution, the virome family and species composition and abundance differed to some extent between the two species, and also by geographic location, at this scale of sampling ( Figure 1, Figure 2, Table 2).
Viewed at the family/order level, the relative virome abundance of Cx. pipiens was dominated by luteo-, orthomyxo-, and the Nam Dinh nidovirus. In comparison, Cx. torrentium was dominated by the Nam Dinh nidovirus in addition to the picorna-, mononega-, and reo-viruses. It should be noted, however, that family-wide comparisons could be skewed by the presence of single highly abundant viruses, as was the case here (particularly Nam Dinh virus that represented 30% of all reads in library 2), such that analyses of relative abundance and diversity are better conducted at the species level. Viewed at the level of species per region, we identified several viruses that were seemingly unique to their respective sampling location ( Figure 3). This suggests that local acquisition, as well as local ecosystem and habitat composition, may be important in shaping virome compositions, although this will need to be confirmed with a larger sample size of mosquitoes. Although we did not find any evidence that sampling year impacted the results, as the majority of viruses were found in libraries covering different sampling years ( Table 2, Table S2) it is likely that detailed longitudinal sampling would provide more information on possible seasonal changes in virome composition.
Direct comparisons between published virome studies are complicated by such factors as differences in sequencing technologies, bioinformatic analyses, criteria for species demarcation, and study focus. Despite these important caveats, it is noteworthy that the number of viruses found in the relatively small sample here is of a similar magnitude and diversity to those found at lower latitudes [12,15]. Hence, the virome composition appears not to follow the same trend as mosquito biodiversity, with fewer species in temperate regions [17,18]. Specifically, 24 different viruses were found in Cx. torrentium, of which six were unique to that species, and 34 viruses were found in Cx. pipiens, of which 16 were unique. Hence, 18 viruses were shared between both mosquito species, 16 of which we tentatively consider to be mosquito associated based on their abundance and phylogenetic position (Figures 3-6, Table 3).
Given their relatively close phylogenetic relationship (Figure 7), and the fact that both mosquito species inhabit the same region, share larval habitat [4], and blood-meal hosts [30], the difference in virome composition between Cx. pipiens and Cx. torrentium is striking. By considering virus abundance and phylogenetic position we suggest that 26 of the viruses discovered were likely mosquito associated (Figures 4-6, Table 3), although we cannot exclude either false-negative or false-positive associations. For example, the divergent Ista virus (Picornaviridae) was found in high abundance and in multiple libraries but did not cluster with any viruses that originated from mosquitoes, although it did group with other arthropods (Figures 4-6). The fact that it did not cluster with other mosquito viruses is perhaps unsurprising as studies from temperate regions are few, and this is the first study investigating the virome of Cx. torrentium. It is clear that many viruses are seemingly ubiquitous in mosquitoes, covering a wide variety of climates and habitats [12,15,31], but whether Ista virus and many other viruses are truly mosquito associated will need to be considered in additional studies. In particular, virus isolation and cell-culture/lineage experiments will be central to determining the host association of the Ista virus and other viruses discovered via meta-transcriptomic studies. It was also noteworthy that no insect-specific flaviviruses were discovered in this study, even though these are relatively commonplace [32] and have previously been found in mosquitoes in Northern Europe [33,34]. Relatedly, although we found no evidence that the virus sequences obtained here are from endogenous viral elements (EVEs) [35,36], we cannot definitively exclude that some of the viruses documented may in fact be derived from other organisms and/or host genomes present within or on the outside of the mosquito. Notably, our study reveals that pathogenic viruses such as SINV can sometimes have similar abundance levels to viruses not associated with human disease (Table 2), in turn suggesting that pathogenic viruses form a natural part in the overall virome composition. The difference in the abundance of specific viruses between Cx. torrentium and Cx. pipiens is interesting. It has previously been shown that Cx. pipiens is commonly infected with Wolbachia, while this bacterium is absent from Cx. torrentium [9]. Wolbachia is well-known for its ability to block virus infection in some mosquito species, although it has mostly been studied in systems with pathogenic viruses such as dengue [37]. Although one potential explanation for the difference in virome composition is differential associations with the intracellular bacteria Wolbachia pipientis, we found no compelling evidence for Wolbachia in any of the Culex samples studied here.
The species separation between Cx. pipiens and Cx. torrentium has long been ignored, largely because of the need for molecular demarcation, so that it has been assumed that most of the biology of the two species is comparable. Our study indicates that these two species have differing virome compositions and also that Culex mosquitoes in northern temperate regions can harbor similar viral diversity as mosquitoes in tropical and sub-tropical regions. Further studies should consider the host  Figure 7. Phylogenetic relationships, based on partial COX1-gene, of Cx. pipiens and Cx. torrentium for all libraries (L1-L12) together with representative publicly available reference sequences (with their associated GenBank accession numbers). Numbers on branches indicate SH support, and only branches with SH support ≥80% are indicated. Branch lengths are scaled according to the number of nucleotide substitutions per site. The tree is midpoint rooted for clarity only.
Notably, our study reveals that pathogenic viruses such as SINV can sometimes have similar abundance levels to viruses not associated with human disease (Table 2), in turn suggesting that pathogenic viruses form a natural part in the overall virome composition. The difference in the abundance of specific viruses between Cx. torrentium and Cx. pipiens is interesting. It has previously been shown that Cx. pipiens is commonly infected with Wolbachia, while this bacterium is absent from Cx. torrentium [9]. Wolbachia is well-known for its ability to block virus infection in some mosquito species, although it has mostly been studied in systems with pathogenic viruses such as dengue [37]. Although one potential explanation for the difference in virome composition is differential associations with the intracellular bacteria Wolbachia pipientis, we found no compelling evidence for Wolbachia in any of the Culex samples studied here.
The species separation between Cx. pipiens and Cx. torrentium has long been ignored, largely because of the need for molecular demarcation, so that it has been assumed that most of the biology of the two species is comparable. Our study indicates that these two species have differing virome compositions and also that Culex mosquitoes in northern temperate regions can harbor similar viral diversity as mosquitoes in tropical and sub-tropical regions. Further studies should consider the host range of these viruses, their potential interactions with pathogenic viruses, and how virome composition is determined by mosquito host structure and feeding preference. Isolating the viruses discovered here will also be a key priority to enable a better understand of mosquito-virus co-evolution. Ultimately, it will be essential to identify the key evolutionary, ecological, and environmental factors that determine virome composition, as well as the impact of virome composition on the mosquito.
Supplementary Materials: The following are available online at http://www.mdpi.com/1999-4915/11/11/1033/s1, Table S1: Information on collection site, year of collection and pool size for all Cx. pipiens and Cx. torrentium libraries; Table S2: Number of reads mapped to each virus, Wolbachia and host genes per library per mosquito species.