Virome of Mosquitoes from Natural Landscapes of the Western Siberia

1. Introduction

Mosquitoes are well known to transmit numerous arboviruses causing viral infections in animals and humans, such as West Nile virus (WNV), Zika, Japanese encephalitis, Chikungunya, and dengue viruses [1,2,3,4]. In recent years, the development of metagenomic approaches has led to the discovery of many novel viruses in invertebrates [5,6,7,8,9,10]. Studying the viromes of different species of mosquitoes has revealed new viruses referred to as Insect-Specific Viruses (ISVs). The viral interference within the ISV group and pathogenic viruses may dramatically change the viral biodiversity in mosquitoes and thereby predetermine the transmission of pathogenic viruses by mosquitoes to animals and humans [11,12,13,14,15,16]. These are viruses belonging to Peribunyaviridae, Flaviviridae, Reoviridae, and Togaviridae families, with all of them being a potential source of viral biodiversity for viruses with dual tropism for invertebrate and vertebrate hosts [7,17,18].

Comparing mosquito viromes from different geographical regions revealed their biodiversity, providing new insights into the phylogeography of mosquito-borne viruses [5,6,19,20]. Generally, such studies are conducted in countries characterized by warm or tropical climates, such as China, Australia, Mozambique, and the USA, where mosquito-borne viral infections are not uncommon. Sindbis, Inco, and West Nile viruses are usually detected and isolated from mosquitoes in the southern regions of Russia [21,22]. No systemic information is available for Western Siberia, where it is also possible for mosquito-borne viruses to circulate. This region has a continental climate with long winters and short summers that may limit the biodiversity of mosquito species and mosquito-associated viruses. Only several studies have reported the detection of WNV markers in birds and human cases of West Nile fever in Western Siberia [22,23].

In this study, we sought to investigate the biodiversity of mosquito-borne viruses from different mosquito species collected in Western Siberia using metagenomic approaches.

2. Materials and Methods

2.1. Mosquito samples

For the study, 3.910 mosquitoes were collected in the Novosibirsk region during the spring-summer period of 2017–2018. The collection sites were selected in typical mosquito habitats in Western Siberia (Figure 1.): deciduous and mixed forests with well-developed grassy cover, deforestations with a natural resumption of hardwoods, and stream banks. The mosquitoes were collected using a light trap [https://survinat.ru/2011/09/metodika-sborov-xraneniya-i-izucheniya-komarov/]. The capture was conducted after sunset. The mosquitoes were transported in a thermal bag, on a damp napkin, at a temperature of 4^°C and stored at minus 18–24^°C. The fragments of the 16S rRNA and COI gene of the mitogenome were sequenced to determine the mosquito species [11]. The pools of 10–40 mosquitoes were formed using the data on mosquito species and the time of collection.

2.2. Sample preparation

All the mosquitoes were washed in 70% ethanol and then twice in water, followed by homogenization to remove potential surface microorganisms. The homogenization of the samples was performed mechanically by grinding in a mortar with 300 µl of sterile saline. The homogenates were centrifuged at 8,000 g for 5 minutes at 4^°C, and the supernatants were used for the analysis. The total RNA was extracted using Extract RNA reagent (Eurogen, Russia) according to the manufacturer’s protocol and purified on Cleanup Mini spin columns (Eurogen, Russia). The pools were then processed by benzonaze [24]. The first chain of cDNA was synthesized using the NEBNext Ultra Direction module. The second cDNA chain was synthesized using the UMI Second Strand Synthesis module (Illumina, Lexogen).

2.3. NGS sequencing and phylogenetic analysis

The dsDNA libraries were prepared and analyzed by NGS on MiSeq using the Illumina technology. Cutadapt (version 1.18) and SAMtools (version 0.1.18) were used to remove Illumina adapters and re-read. The contigs were assembled de novo using the MIRA assembly (version 4.9.6). The experimentally determined sequences were deposited in GenBank. The phylogenetic analysis was performed using RNA-dependent RNA polymerase (RdRp) sequences from GenBank with amino acid identity > 20%. The sequences were aligned, and phylogenetic trees were built in Vector NTI Advance 11, MEGA 7/10 (PSU, USA), and Lasergen 7 (Invitrogen). The resulting viral sequences and sequence read archive (SRA) were deposited in GenBank. The phylogenetic trees were calculated by the maximum likelihood method using 500 replicates for bootstraps values.

3. Results

3.1. Mosquito species

A total of 3,910 mosquitoes were collected in the Novosibirsk suburbs and Novosibirsk region rural district in 2017–2018. The pools of 10–40 mosquitoes were formed to identify the mosquito species for every collection point. The fragments of the 16S rRNA and COI gene of the mitogenome were sequenced to determine the species of mosquitoes (Figure 2). Aedes caspius (Pallas, 1771) and Ae. mariae (Sergent and Sergent, 1903) were found to be the most abundant species, making up to 41.2% of the total. Anopheles messeae (Falleroni, 1926) accounted for 34.6%, Culex pipiens (Linnaeus, 1758) for 8.3%, Coquillettidia richardii (Ficalbi, 1889) for 5.5%, C. modestus (Ficalbi, 1889) for 4.8%, An. maculipennis (Meigen, 1818), and An. sinensis (Wiedemann, 1828) for 2.8%.

3.2. NGS sequencing

One hundred forty-four putative viral sequences were selected in the first step, and of these, 30 sequences with lengths greater than 1239 bp were chosen. Eight sequences with a level of identity > 80% aa) have been previously described as mosquito-borne viruses (Table 1). These are Partitivirus-like 1 (dsRNA, Partitiviridae), Hammarskog tombus-like virus (ssRNA (+), Tombusviridae), Hammarskog picorna-like virus (ssRNA (+), Picornavirales, unclassified), Lymantria dispar iflavirus 1 (ssRNA (+) viruses, Iflavirus), Wenzhou noda-like virus 6 (ssRNA (+) viruses, unclassified), Mayapan virus (ssRNA (+), 2 segments, Sanxia permutotetra-like virus 1 (ssRNA (+) viruses, unclassified), Chaq virus-like 1 (RNA viruses, unclassified). Other 22 viral sequences are presented as putative mosquito-borne sequences, with a level of identity less than 79% for viral prototype sequences.

The proportions for classified and unclassified viral reads are presented for 13 mosquito pools (Figure 3). The prevalence of unclassified and classified picornaviruses was detected practically in all studied pools. Unclassified viral sequences were also analyzed, ranging from 0.13 to 34.27% for the different pools.

The phylogenetic analysis results for the sequences obtained from mosquitoes are presented in Figure 4. This phylogeny is based on an amino acid sequence of RdRp, with these data confirming the virome biodiversity of mosquito viruses in nature.

3.2.1. Tymovirales (positive ssRNA)

The complete viral genome tymovirus-like sequence with a 60% identity (according to amino acid sequence) with the previously described Insect-associated tymovirus 1 in Mexico (MN203215) was detected in the Cq. richardii mosquitoes. This virus has been designated as Inya insect-associated virus 1 (MW251313, MW251314). The genomic positive ssRNA of the Inya insect-associated virus was identified to comprise 6526 bp and three ORF-encoding proteins (Figure S1). ORF MP of Inya insect-associated virus 1 is RdRp, and this ORF contains a highly conservative "tymobox" near the 5’-end [25]. The tymobox sequence has 16 nucleotides that are likely part of the subgenomic promoter for the third ORF encoding the coat protein (CP). Previously, the Tymovirales were well-known as plant viruses [26]. Inya insect-associated virus 1 can be presumably identified taxonomically by the order Tymovirales, unclassified Tymovirales.

3.2.2. Partitiviridae (dsRNA)

Partitivirus-like 1 was detected in Cq. richardii with an identity level of 89% with an isolate from Anopheles gambiae collected in Liberia (KX148575). Another novel putative partitivirus was detected in the pool of Cq. richardii mosquitoes (Figure S2) and was designated as Krahall insect-associated virus 1 and 2 with an identity level of 60–63% with previously described Atrato Partiti-like virus 2 isolated earlier from Anopheles darlingi in Colombia (MN661058). In addition, seven suspected partitiviruses were found in the pool of Cq. richardii mosquitoes. These are: novel insect Talaya 1 and 2 viruses (MW251327–MW251330), with a 68–71% aa identity with previously described Partitivirus-like 1 (Liberia, KX148575); Tarbrook virus (MW251325) with 79% homology with the previously described Sonnbo virus in Sweden (MK440649); Zeyabrook partiti-like_viruses 1 and 2 having similarity of 38–42% with Beihai partiti-like virus 2 (NC_032500) from China. All the prototype sequences were isolated earlier from invertebrates (mollusks, octopuses, mosquitoes, and odonatos). These partitiviruses were preliminarily taxonomically identified as the family Partitiviridae, unclassified Partitiviridae.

3.2.3. Totiviridae (monopartite dsRNA)

We have found a novel putative totivirus designated as Zyryana toti-like virus 2 with the prototype Fitzroy Crossing toti-like virus 1 isolated earlier from Culex annulirostris in Australia (MT498830) (Table 1). Zyryana toti-like virus 2 with a 51% identity with prototype sequence was detected in the pools of An. Messeae and Cq. richardii mosquitoes. The length of the nucleotide sequence of the Zyryana toti-like virus 2 was over 5863 bp, and the ORFs encode two proteins 1112 aa and 807 aa of CP and RdRp. The putative genome organization schemes for these totiviruses are presented in Figure S3.

3.2.4. Tombusviridae (positive ssRNA).

Tombusviridae (Tolivirales, Tombusviridae) are single-stranded RNA (+) genomes between 3.7 and 4.8 kb in length, currently regarded as plant viruses with a relatively limited host [27]. We presented the complete polyprotein (4317 bp) for Hammarskog tombus-like virus with 90% identity with a similar virus detected in Sweden (MN513379) and isolated from Cq. Richiardii in 2017. The 4166 bp partial polyprotein contains three ORF-encoded hypothetical polypeptides: 397 aa, 482 aa, and 409 aa that differ for Hubei tombus-like 20 (Figure S4). In addition, we found a novel Oyosh tombus-like virus, with 62% identity with Hubei tombus-like virus 13 (NC033017) isolated from house centipedes in China. Four polypeptides are encoded by a prototype genome (5904 bp). The RdRp for Tombusviridae was translated using a potential alternative mechanism to suppress the stop-codon reading mechanism with the formation of a full-size protein with an elongated C-end of ORF1 [6].

3.2.5. Picornavirales (positive ssRNA)

Most members of the order Picornavirales have a single molecule of positive sense RNA ranging in length between 7,000 and 12,500 nt. The viral RNA is infectious and serves as a template for replication and mRNA [28]. Six different picorna-like viruses were identified as mosquito-associated viruses in Western Siberia (Table 1). We have assembled a complete genome for the Hammarskog picorna-like virus (11507 bp) from Cq. richardii mosquito that has five OFR encoding 175 aa, 156 aa, 121aa, 376 aa, and 2424 aa polypeptides with 98% aa identity with the previously described Hammarskog picorna-like virus (MN513381) isolated from Cq. richiardii in Sweden (Figure S5). In addition, other novel picorno-like viruses were found in Cq. richardii mosquitoes collected in the Novosibirsk region. These are Miltyush picorna-like viruses 1 and 2, Isses picorna-like virus, Ichacreek insect virus, and polycipiviridae associated with Ora rivulet insects.

Miltyush picorna-like virus was found to have only 36% identity with the previously detected Halhan virus 3 from Haliotis discus hannai in Korea (NC040628). Isses picorna -like virus 1 was found to have 64% identity with previously discovered Washington bat picornavirus in the USA (KX580885). Ichacreek insect virus 3 was identified to have a 44% level of identity with the previously detected Solenopsis invicta virus 3 (GU017972) from Solenopsis invicta in Argentina. Ora rivulet insect-associated polycipiviridae was identified to have a 34% identity level with the previously discovered Polycipivirida sp. isolated from Pteropus lylei in Cambodia (MK161350). The preliminary taxonomic identification for these six viruses is Picornavirales, unclassified Picornavirales.

3.2.6. Iflaviridae (positive ssRNA)

The order Picornavirales also includes some iflaviruses that were found in Cq. richardii mosquito in this study. Lymantria dispar iflavirus 1 was detected with its sequence having 99% identity with already known viral isolates in USA and Russia (KJ629170, MN938851). The alignment and phylogenetic analysis revealed a high sequence identity with the representatives of Iflavirus, the family Iflaviridae (data not shown).

3.2.7. Nodaviridae (Bi-partite positive-sense, ssRNA)

We have assembled practically the whole genome for Wenzhou noda-like virus 6 from Cq. richardii mosquito and Mayapan virus with 87 and 90% identity levels, respectively. Previously, the sequence of the Wenzhou noda-like virus 6 was identified in Channeled applesnail in China (KX883260) and Mayapan virus (MH719096) isolated from the Psorophora ferox mosquito in Nexico (Figure S6). Other novel nodaviruses were also found in Cq. richardii mosquitoes collected in the Novosibirsk region. These are Mayzas noda-like virus RNA 2 with a prototype Mayapan virus RNA2 segment (MH719097) with a 43% identity with that isolated from Psorophora ferox mosquito in Mexico and the insect-associated Uzakla virus with a prototype Mosinovirus (KJ632942) with a 36% identity and isolated from Culicidae spp. in Cote d’Ivoire.

3.2.8. Permutotetraviridae (dsRNA)

The genomic RNA for permutotetraviruses is 4,582 bp long and encodes three ORFs overlapping in a short region (Figure S7). The longest ORF (1028 aa) encoding RdRp overlaps with 106 nucleotides with a small ORF (199 aa), presumably encoding the capsid protein. Like all permutotetraviruses, the sequences from Cq. richardii mosquito pools showed the presence of the virus with a 93% identity with previously detected Sanxia permutotetra-like virus 1 in water striders in China (KX883450). The Uzakla mosquito-associated permutotetra-like virus with a 53% identity to Vespa velutina permutotetra-like virus 2 in France (MN5650551, MN565052) was early described as unclassified Permutotetraviridae.

3.2.9. Other viruses and viral sequences

Two variants of Chaq virus-like 1, Tartas insect associate virus, ZeyaBrook chaq-like virus 2, Kamenka insect-associated virus, and Uzakla insect virus were detected in Cq. richardii. The Chaq virus-like 1 has 82% identity with an earlier described unclassified sequence from Anopheles gambiae in Liberia (KX148554). The ZeyaBrook chaq-like virus 2, Kamenka insect-associated virus, and Uzakla insect virus have a 33–47% identity with the previously unclassified putative viral sequences (KX148556, KX883594, and NC032218) isolated from invertebrates. Only the Tartas insect associate virus may be classified as unclassified Solemoviridae with prototype Atrato Sobemo-like virus 6 (MN661101) with a 61% identity detected in Wyeomyia spp. mosquitoes in Colombia. The Solemoviridae have a relatively small (4–4.6 kb) positive-sense, single-stranded, monopartite RNA genome with 4–5 ORFs, and they are usually associated with plant viruses.

4. Discussion

The application of the metagenomic approach offers novel opportunities for virome analysis [5,7,20]. This approach has provided new insights into the evolution of viruses of clinical importance and has allowed new viruses to be discovered from different viral families such as Peribunyaviridae [5], Rhabdoviridae [29], Orthomyxoviridae [7,30], Flaviviridae [31] and Reoviridae [32], as well as unclassified Chuvirus [7], and Negevirus [33]. Recent metagenomic studies have also confirmed the presence of dengue virus, Zika virus, and Japanese encephalitis virus in mosquitoes in China [34,35].

Numerous genetically diverse viruses have also been detected by NGS sequencing in plants, invertebrates and vertebrates in tropical countries [7,8]. Phylogenetic analysis has demonstrated that it is possible for all host species and viruses to co-evolve by changing hosts. Mosquitoes are among the most common and important viral vectors of the Zika, dengue, yellow fever, and West Nile viruses that are associated with unprecedented global outbreaks of these infectious in tropical countries [36,37]. In addition, mosquitoes are also known to carry insect-specific viruses. Although not directly affecting humans and animals, these viruses can modulate the transmission of pathogenic viruses to vertebrates [38,39]. The growth of tourism and trade has also led to an intensive exchange of viral pathogens and their vectors in different geographic regions. Together with the rapid growth of large cities in tropical countries, these are the basis for outbreaks or/and epidemics for mosquito-borne infections among animals and humans, with the environment to maintain the transmission of zoonotic infectious [40]. In addition, viruses have extraordinary evolutionary potential to generate new pathogenic isolates that can cause severe diseases in humans and/or animals.

The south of the Western Siberian Plain is characterized by a continental climate, with short warm summers and long winters, uniform humidity, and rather abrupt changes in all-weather components over relatively short periods of time [41]. This region has experienced characteristic negative mean annual temperatures during the last century, with the maximum variations of the mean annual temperature being 3.6 °C over the observation period. The activity season for different species of mosquitoes begins when the ambient temperature rises above 0 °C (early May) and ends in late August or early September, depending on the year. The maximum duration of the mosquito activity period is approximately four to five months. Seventeen species of mosquitoes were earlier found in the forest-steppe and steppe zones of the region [42]. The mosquito species composition from different foci can drastically vary. For example, the Cq. richardii concentrations can vary from 1.7 to 99.5%, with this species usually dominating in the main forest-steppe and steppe landscapes of the rural part of the Novosibirsk region.

In this study, we used a metagenomic sequencing method to identify the viromes in seven mosquito species collected in the vicinity of Novosibirsk. The metagenomic approach was used to identify the viral diversity in randomly collected mosquitoes. We have identified 30 coding complete viral genomes and viral-like partial sequences of capsid proteins and/or RdRp from mosquitoes (Table 1). These sequences were classified as putative members of orders: Tymovirales and Picornavirales, families: Partitiviridae, Totiviridae, Tombusviridae, Iflaviridae, Nodaviridae, Permutotetraviridae, Solemoviridae, and four unclassified RNA-viruses. The previously described Partitivirus-like 1, Hammarskog tombus-like virus, Hammarskog picorna-like virus, Lymantria dispar iflavirus 1, Wenzhou noda-like virus 6, Mayapan virus, Sanxia permutotetra-like virus 1, and Chaq virus-like 1 were identified as practically complete genomes with an 82–99% level of identity in Cq. richardii mosquito. These viruses were earlier found in Liberia (West Africa), Sweden (North Europe), the USA (America), China (Asia), and Mexico (Central America). These findings allow us to hypothesize that these viruses may be widely distributed on a global scale.

Some novel putative viruses and viral sequences have prototype viral sequences with 31% to 79% identity levels, with these prototypes also found in invertebrates from almost all continents. Some of them are associated with different species of mosquitoes. In our study, the main parts of novel viruses were associated with Cq. richardii mosquito, with this species widespread in the south of Western Siberia [42]. The role of this mosquito species in the spread of human viral infections in Siberia has not been studied virtually, suggesting that our knowledge concerning mosquito-associated viruses in North Eurasia is very limited and requires further study.

5. Conclusions

We have identified novel and known viral genomes and viral-like partial sequences in mosquitoes collected in the Novosibirsk region of Western Siberia. They were classified as novel putative viruses using the bioinformatics analysis of partial sequences of capsid proteins and RdRp or whole polyproteins (genomes) within the orders Tymovirales and Picornavirales, the families Partitiviridae, Totiviridae, Tombusviridae, Iflaviridae, Nodaviridae, Permutotetraviridae, and Solemoviridae, and four unclassified RNA-viruses. We believe that the virus identification will enhance our understanding of the transmission of RNA viruses by mosquitoes in North Asia. We hope that the discovery and observation of these mosquito-borne viruses can help prevent future outbreaks of viral infections in the region under study.