Mosquito-borne viruses, insect-specific flaviviruses (family Flaviviridae, genus Flavivirus), Banna virus (family Reoviridae, genus Seadornavirus), Bogor virus (unassigned member of family Permutotetraviridae), and alphamesoniviruses 2 and 3 (family Mesoniviridae, genus Alphamesonivirus) isolated from Indonesian mosquitoes.

Mosquitoes transmit many kinds of arboviruses (arthropod-borne viruses), and numerous arboviral diseases have become serious problems in Indonesia. In this study, we conducted surveillance of mosquito-borne viruses at several sites in Indonesia during 2016-2018 for risk assessment of arbovirus infection and analysis of virus biodiversity in mosquito populations. We collected 10,015 mosquitoes comprising at least 11 species from 4 genera. Major collected mosquito species were Culex quinquefasciatus, Aedes albopictus, Culex tritaeniorhynchus, Aedes aegypti, and Armigeres subalbatus. The collected mosquitoes were divided into 285 pools and used for virus isolation using two mammalian cell lines, Vero and BHK-21, and one mosquito cell line, C6/36. Seventy-two pools showed clear cytopathic effects only in C6/36 cells. Using RT-PCR and next-generation sequencing approaches, these isolates were identified as insect flaviviruses (family Flaviviridae, genus Flavivirus), Banna virus (family Reoviridae, genus Seadornavirus), new permutotetravirus (designed as Bogor virus) (family Permutotetraviridae, genus Alphapermutotetravirus), and alphamesoniviruses 2 and 3 (family Mesoniviridae, genus Alphamesonivirus). We believed that this large surveillance of mosquitoes and mosquito-borne viruses provides basic information for the prevention and control of emerging and re-emerging arboviral diseases.

Mosquitoes transmit many kinds of arboviruses belonging to families such as Flaviviridae, Togaviridae, Phenuiviridae, and Rhabdoviridae [2]. These mosquito-borne arboviruses cause endemics or epidemics in tropical/subtropical regions worldwide, depending on the distribution and activities of vector mosquitoes [56]. In addition, a recent study showed that mosquitoes also harbor a diverse range of RNA viruses that do not seem to infect human and other vertebrates. These viruses are known as insect-specific viruses (ISVs) and can infect and replicate only in invertebrates [45]. ISVs have been identified in several viral taxa such as Flaviviridae, Phenuiviridae, Mesoniviridae, negeviruses, Reoviridae, Nodaviridae, Rhabdoviridae, and Togaviridae [6, 44, 45]. Because some reports have revealed that several ISVs have been found in immature stages of mosquitoes [4, 36, 50], they have been probably maintained by vertical transmission from infected female mosquitoes [3]. The study of these various ISVs in vector mosquitoes would contribute to the understanding of virus biodiversity and the complex vector-borne viral disease system in nature because ISVs have the potential to disrupt arbovirus transmission or reduce vector abundance [6].

Although there are many factors for global emergence and re-emergence of arboviral diseases, one of them is change in the environment surrounding human societies, such as urbanization of rural areas by human migration, global tourism, and expansion of agriculture areas. In addition, trade and transportation have contributed to the introduction of both arboviruses and vector arthropod species. For instance, dengue virus (DENV; family Flaviviridae, genus Flavivirus) is one of the most common mosquito-borne arboviruses and more than 2.5 billion people are at risk of this viral disease worldwide [18]. DENV has been transmitted by Aedes mosquito species, such as Aedes aegypti (Linnaeus) and Aedes albopictus (Skuse), and distributions of these mosquito species have been strongly influenced by urbanization factors due to their preference of artificial containers as its larval habitat and internationalization [18]. Climate change is also associated with emergence and re-emergence of arbovirus diseases because the distribution of arboviruses is dependent on their arthropod vectors [40].

In Indonesia, many arboviral diseases are prevalent throughout the year. Recently, 59,047 cases of dengue fever, including 444 deaths, were reported in Indonesia during 2016–2017 [38]. Moreover, 326 and 1,702 cases of Japanese encephalitis virus (JEV; family Flaviviridae, genus Flavivirus) and Chikungunya virus (CHIKV; family Togaviridae, genus Alphavirus), respectively, were reported during 2016 in Indonesia [38]. Moreover, Zika virus (family Flaviviridae, genus Flavivirus) infection was first reported in Jakarta in 2012 [31]. The number of patients with arbovirus infections has tended to increase annually, as well as the number of mosquito vectors and distribution of arboviruses in Indonesia [54].

Entomological approaches may aid in risk assessment of arbovirus infections. In this study, we conducted surveillance of mosquito species and mosquito-borne viruses in Indonesia during 2016–2018.

MATERIALS AND METHODS

Mosquito collection

Mosquito collection was conducted from seven provinces of four islands in the eastern part of Indonesia (Bengkulu province in Sumatra; Banten, Jakarta, West Java, and Central Java provinces in Java; Bali province in Bali; and Central Kalimantan province in Borneo) during 2016–2018 (Fig. 1). Mosquitoes were collected in animal houses, such as swinery (Bengkulu and Denpasar), cowsheds (Bogor), and non-human primate sanctuaries (Palangka Raya and Bogor), and both indoor and outdoor human residential areas (Bogor, Tangerang, Bengkulu, Solo, Bekasi, and Pasar Minggu). Adult mosquitoes were collected in residential areas at daytime and around animal houses after sunset by using sweep nets and aspirators. Blood-feeding mosquitoes were kept in cages for a few days to let them digest the blood in their guts. The collected mosquitoes were identified based on morphology [48]. After identification, mosquitoes were sorted into pools with a maximum of 96 adults per pool and kept at −80°C for virus isolation.

Fig. 1.

Mosquito collection sites in this study. Black points show Bengkulu province (1) in Sumatra; Banten (2), Jakarta (3), West Java (4), and Central Java (5) provinces in Java; Bali province (6) in Bali; and Central Kalimantan province (7) in Borneo.

Cell culture

Two mammalian cell lines, Vero [Japanese Collection of Research Bioresources (JCRB), Osaka, Japan] and BHK-21 (JCRB), and one mosquito cell line, C6/36 (Health Science Research Resources Bank, Osaka, Japan), were used for virus isolation. Vero and BHK-21 cells were cultured in Dulbecco’s modified eagle’s medium (DMEM; Sigma-Aldrich, St. Louis, MO, USA) containing heat-inactivated 5% fetal bovine serum (FBS; Hyclone Laboratories, South Logan, UT, USA) and 1% penicillin and streptomycin (Pen-Strep; Life Technologies, Carlsbad, CA, USA), and C6/36 cells were cultured in Eagle’s minimum essential medium (Sigma-Aldrich) containing heat-inactivated 10% FBS and 1% Pen-Strep. Cell lines were maintained at 37°C for mammalian cells and 28°C for mosquito cells under 5% CO2.

Virus isolation

Virus isolation was performed as previously described [28, 29]. Each pool of mosquitoes was homogenized in 500 µl of ice-cold DMEM containing 2% heat-inactivated FBS and 1% Pen-Strep using the Multi-beads shocker (Yasui Kikai Inc., Osaka, Japan). Homogenates were clarified by centrifugation (21,500 × g for 5 min at 4°C), and supernatants were passed through sterile 0.45-µm filters (Corning Inc., Corning, NY, USA). Filtrates were inoculated onto a monolayer of Vero, BHK-21, or C6/36 cells in 24-well culture plates (Sumitomo Bakelite Inc., Tokyo, Japan), and plates were incubated for 2 hr to allow virus adsorption. After adding fresh medium, cell cultures were incubated at cell culture conditions. Cytopathic effects (CPEs) were observed daily, and supernatants were collected from cells with CPEs. If no CPE was observed, cells were blind-passaged, and after five blind passages, culture supernatants were collected even if no CPEs were observed. These culture supernatants were stored at −80°C until use.

Detection of viral genome

RNA was extracted from culture supernatants using the QIAamp Viral RNA Mini Kit (QIAGEN, Inc., Valencia, CA, USA). RNA was extracted from each culture supernatant from cells with CPEs individually. By contrast, for cells without any CPEs, culture supernatants from 5–8 wells of Vero and BHK-21 cells or 3–4 wells of C6/36 cells were pooled for RNA extraction, which was subsequently used in RT-PCR analyses. If RNA from mixed pools yielded positive RT-PCR results, then RNA was again extracted from each supernatant in the pool and analyzed individually. RT-PCR was conducted using the QIAGEN OneStep RT-PCR kit (QIAGEN) to detect viral genes. Four universal primer sets were used for detection of flaviviruses [26, 51], alphaviruses [12], phleboviruses, and rhabdoviruses [30] (Table 1). RT-PCR for each primer sets were performed at 50°C for 30 min (reverse transcription), 95°C for 15 min, followed by 35 cycles of 94°C for 30 sec, 55°C for 30 sec, and 72°C for 1 min (amplification), and final extension at 72°C for 10 min. PCR products were separated by electrophoresis on 2% agarose gels, and intensive bands were purified by using the MinElute Gel Extraction Kit (QIAGEN). Purified PCR products were directly cycle-sequenced using the ABI PRISM BigDye Terminator Cycle Sequencing Kit v.3.1 (Thermo Fisher Scientific, Waltham, MA, USA) and ABI PRISM 3130 Genetic Analyzer (Thermo Fisher Scientific). Minimum infection rate (MIR) was calculated as the ratio of the number of positive pools to the total number of mosquitoes tested. For samples from cells that showed clear CPEs but yielded negative RT-PCR results for viral genes, high-throughput analyses were applied to identify viral genomes as described below.

Table 1.

List of primer sets used in this study

Primer name	Direction	Sequence	Size (bp)	Target virus	Target gene	Reference
MAMD	Sense	5′-AACATGATGGGRAARAGRGARAA-3′	260	Flavivirus	NS5	[51]
cFD2	Antisense	5′-GTGTCCCAGCCGGCGGTGTCATCAGC-3′				[26]
VIR2052F	Sense	5′-TGGCGCTATGATGAAATCTGGAATGTT-3′	144	Alphavirus	nsP4	[12]
VIR2052R	Antisense	5′-TACGATGTTGTCGTCGCCGATGAA-3′
MBPL3100F	Sense	5′-AGTCTCYTCTGCCATYTC-3′	890	Phlebovirus	L	This study
MBPL3287R	Antisense	5′-AGGATCTRGARGGGAACTTRT-3′
RHNB1520F	Sense	5′-ACIAAIAARTWIATGATGATGAA-3′	188	Rhabdovirus	N	[30]
RHNB2134R	Antisense	5′-TGIARDATICCYTGCATCAT-3′
BANV-F	Sense	5′-AGATCCTAACTGTGACCCAATGTT-3′	770	Banna virus	VP1	This study
BANV-R	Antisense	5′-TGTAACTTCTAACAAATCCGCAAA-3′
BGV-F	Sense	5′-GTAGACGAATGCATGTTCGATAAG-3′	441	Bogor virus	putative RdRp	This study
BGV-R	Antisense	5′-CCGTCTAACTGTGTGGATAACAAG-3′
AMSV-F	Sense	5′-TATGGCAAACGACGTATAGCAG-3′	371	Mesonivirus	putative RdRp	This study
AMSV-R	Antisense	5′-AAGCATARAYTGGTTGTGACG-3′

High-throughput analyses

Three selected culture supernatants of C6/36 cells showing different forms of CPEs (pool nos. 2, 34, and 69) were subjected to next-generation sequencing (NGS) analyses (Supplementary Fig. 1). Approximately 50 ml of each culture supernatant were concentrated to approximately 0.5 ml by using Amicon Ultra 50 kDa (Merck Millipore, Billerica, MA, USA). Concentrated supernatants were ultracentrifuged at 60,000 × g for 3 hr at 4°C through a 25% sucrose cushion. Then, supernatants were carefully removed and pellets were resuspended in 200 µl of sterile PBS on ice. Total RNAs were extracted from samples by using the QIAamp Viral RNA Mini Kit (QIAGEN) according to the manufacturer’s instructions. We applied three methods for NGS analyses: Ion Torrent PGM system (Life Technologies), Illumina MiSeq system (Illumina, San Diego, CA, USA), and MinION nanopore sequencing device (Oxford Nanopore Technologies, Oxford, England, UK).

For the Ion Torrent PGM analysis, we used RNA obtained from the culture supernatant of C6/36 cells inoculated with mosquito pool no. 2 (Ae. albopictus collected in Bogor, 2016) as previously described [24, 43]. Briefly, RNA was reverse-transcribed into cDNA, which was then amplified with the SeqPlex RNA Amplification Kit (Sigma-Aldrich). The product was end-repaired and ligated with a barcoded adapter using the Ion Xpress Barcode Adapter Kit (Life Technologies) by T4 ligase (New England Biolabs, Ipswich, MA, USA). Barcoded fragments were subjected to emulsion PCR using the Ion OneTouch 2 system with the Ion PGM Template OT2 200 Kit (Life Technologies). Amplified library beads were subjected to the Ion Torrent PGM sequencer (Life Technologies) with the Ion 318 Chip (Life Technologies). Reads obtained by the Ion Torrent PGM were assembled with CLC Genomics Workbench 7.0 software (CLC bio/QIAGEN, Aarhus, Denmark). BWA-MEM [32] was used to map the NGS reads on reference of the 12 genome segments (NC_004198, NC_004200, NC_004201-4, NC_004211, NC_004217-21). SAMtools [33] was used to calculate the coverage.

For the Illumina MiSeq analysis, we used RNA obtained from culture supernatant of C6/36 cells inoculated with mosquito pool no. 34 (Ae. albopictus collected in Bogor, 2017). The DNA library was prepared by using the NEB Next Ultra RNA Library Prep Kit for Illumina (New England Biolabs). Adaptor ligation was conducted by TruSeq 96 dual-index adapters (Illumina) instead of adapters supplied in the kit. Sequencing of paired ends was performed using MiSeq Reagent Kit v2 (300 cycles) (Illumina). Reads were assembled and large contigs were constructed with CLC Genomics Workbench 7.0 software (CLC Bio/QIAGEN).

For nanopore MinION sequencing, we used RNA obtained from the culture supernatant of C6/36 cells inoculated with mosquito pool no. 69 (Cx. tritaeniorhynchus collected in Bogor, 2017). RNA was reverse-transcribed into cDNA, which was amplified with the SeqPlex RNA Amplification Kit (Sigma-Aldrich). After purification using the QIAquick PCR Purification Kit (QIAGEN), synthesized DNAs were directly sequenced by MinION (Oxford Nanopore Technologies) with the Rapid Sequencing Kit (Oxford Nanopore Technologies) according to the manufacturer’s instructions. Reads obtained by MinION were subjected to homology search (BLASTX analysis) using the National Center for Biotechnology Information (NCBI) database.

Phylogenetic analyses

To analyze phylogenetic relationships of isolated viruses, available sequences of related virus species or strains were obtained from the NCBI database. Analysis of the aligned matrix data and the construction of phylogenetic trees were performed using MEGA7 [25] under maximum likelihood algorithms. The best-fit model of nucleotide sequences was found by the model selection analysis by jModelTest 2 [9]. The statistical significance of resulting trees was evaluated by a bootstrap test with 100 replications [11]. Sequence data obtained in this study were deposited in the DNA Data Bank of Japan under accession numbers LC536082–LC536146.

RESULTS

Mosquito collections

A total of 10,015 mosquitoes comprising at least 11 species of 4 genera were collected at several sites in Indonesia in this study (Table 2). The major collected mosquito species were Culex quinquefasciatus (46.4%), Aedes albopictus (20.4%), Culex tritaeniorhynchus (15.8%), Aedes aegypti (10.0%), and Armigeres subalbatus (3.6%). Other minor species, such as Culex hutchinsoni, Culex gelidus, Culex vishnui, Culex fuscocephala, Anopheles vagus, and Anopheles barbirostris, were also collected (Table 2).

Table 2.

Number of mosquitoes used for virus isolation in this study

We compared the species abundance of collected mosquitoes between residential areas and animal houses (Table 2). Although female mosquitoes of Ae. aegypti and Ae. albopictus were collected in some animal houses, they were mainly captured in residential areas, whereas male mosquitoes of Ae. aegypti and Ae. albopictus were captured only in residential areas. For Ar. subalbatus, only female mosquitoes were collected in this study, and they also tended to be primarily captured in residential areas than in animal houses. Cx. quinquefasciatus is one of the common mosquito species in Asian countries, and female mosquitoes of this species were captured in both residential areas and animal houses. By contrast, male mosquitos of Cx. quinquefasciatus were primarily captured in residential areas compared with animal houses. More than 1,500 Cx. tritaeniorhynchus mosquitoes, known as the main vector for JEV, were collected in this study, and approximately 80% of them were collected in one event (Table 2; Bogor, July 2018). Furthermore, female mosquitoes were collected at about 2-fold higher rates than male mosquitoes. Other mosquito species, except for An. vagus, were collected only in residential areas.

The collected mosquitoes were divided into 285 pools and used for virus isolation (Table 3). Seventy-two pools showed clear CPEs in C6/36 cells (25.3%), but not in mammalian Vero and BHK-21 cells. The minimum infection rate (MIR) in C6/36 cells was 0.72 (Table 3). The main manifestation of CPEs in C6/36 cells was cell degeneration, which ranged from swelling, shrinking, and rounded cells to clustering, syncytium formation, and destruction of monolayer cells (Supplementary Fig. 1, represented by pool nos. 2 and 34). Manifestations of CPEs included slower growth rate and increasing gaps of cells. Some CPEs in C6/36 cells formed grape-like clusters 4–6 days postinfection or cells were smaller 3–6 days postinfection (Supplementary Fig. 1, represented by pool no. 69).

Table 3.

Summary of virus isolation (Cytopathic effect (CPE) and RT-PCR) on C6/36 cells from collected mosquitoes

Mosquito species	Sexa)	No of pools tested	CPEs on C6/36	Positive number of RT-PCR [sample no.]
				Flavivirus	Alphavirus	Phlebovirus	Rhabdovirus	Banna virus	Bogor virus	Mesonivirus
Aedes aegypti	F	30	11	4 CFAV [3]	0	0	0	2 [5]	5 [3, 5, 10]	0
Aedes aegypti	M	6	2	2 CFAV [3]	0	0	0	0	0	0
Aedes albopictus	F	33	13	1b) AEFV [3]	0	0	0	4 [1, 5]	7 [3, 5, 10, 18]	2 [8]
Aedes albopictus	M	19	9	0	0	0	0	0	9 [3, 5, 10]	0
Aedes sp.	M	1	0	0	0	0	0	0	0	0
Culex quinquefasciatus	F	99	9	1 CxFV [12]	0	0	0	0	6 [3, 8]	2 [5, 8]
Culex quinquefasciatus	M	10	1	0	0	0	0	0	0	1 [8]
Culex tritaeniorhynchus	F	32	19	0	0	0	0	2 [5]	0	17 [5, 8, 12, 13]
Culex tritaeniorhynchus	M	14	0	0	0	0	0	0	0	0
Culex hutchinsoni	F	2	0	0	0	0	0	0	0	0
Culex vishnui	F	1	0	0	0	0	0	0	0	0
Culex gelidus	F	1	0	0	0	0	0	0	0	0
Culex fuscocephala	F	1	1	0	0	0	0	0	1 [12]	0
Culex sp.	F	1	0	0	0	0	0	0	0	0
Armigeres subalbatus	F	26	7	0	0	0	0	0	6 [8, 10]	1 [5]
Anopheles vagus	F	4	0	0	0	0	0	0	0	0
Anopheles barbirostris	F	4	0	0	0	0	0	0	0	0
Anopheles sp.	F	1	0	0	0	0	0	0	0	0
Total		285	72	8	0	0	0	8	34	23

a) M, male; F, female. b) CPE was not observed.

Detection of viral genomes

RNA was extracted from 72 culture supernatants that showed CPEs in C6/36 cells individually, and the other 3–4 supernatants (213 pools, without CPEs) were mixed into one tube for RNA extraction (total of 132 pools, 72 plus 60 samples for C6/36). Because culture supernatants of Vero and BHK-21 cells did not show any CPEs in this study, 5–8 supernatants were mixed in one tube for RNA extraction (40 samples for Vero and 43 samples for BHK-21). These RNAs (total of 215 RNA samples) were analyzed by RT-PCR using five universal primer sets to detect viruses (Table 1). Only eight samples yielded positive results, which were shown only for the flavivirus primer set; all other samples were negative for all RT-PCR analyses (Table 3). Among the eight flavivirus-positive samples, seven samples were obtained from samples that showed CPEs in C6/36 cells and the other sample was obtained from the mixed Ae. albopictus pool that did not show clear CPEs (Table 3). RNA was then extracted from individual supernatants in the mixed pool and RT-PCR was repeated to determine the virus-positive pool.

Sequencing of these eight RT-PCR products revealed that they showed high homology with sequences of insect-specific flaviviruses in the NCBI database. Six Ae. aegypti (both male and female) showed 92.4–98.0% nucleotide identity to that of Cell fusing agent virus (CFAV) strain Culebra (GenBank accession no. AH015271.2) [8], one Ae. albopictus (female) had 92.4% nucleotide identity to that of Aedes flavivirus (AEFV) strain 21TN (GenBank accession no. KM871198.1) [17], and one Cx. quinquefasciatus (female) demonstrated 96.2% nucleotide identity to that of Culex flavivirus (CxFV) strain Surabaya-2 (GenBank accession no. AB639348.1) [22]. Among these eight RT-PCR-positive samples, CFAV and CxFV showed CPEs in C6/36 cells, but AEFV did not. All insect-specific flaviviruses were obtained from mosquito pools collected in Bogor, West Java province (site no. 4 in Fig. 1 and Table 3).

Only 7 of 72 samples were determined as insect flaviviruses (CFAV and CxFV), and the other 65 CPE-positive samples were not, by RT-PCR analysis with virus-specific primers. To identify CPE-causative agents that were negative for RT-PCR using consensus primers, we utilized the NGS approach.

Detection and isolation of Banna virus

First, we selected one pool (pool no. 2 from Ae. albopictus collected in Bogor, 2016) that showed CPEs in C6/36 cells (Supplementary Fig. 1) for NGS using Ion Torrent (Thermo Fisher Scientific), which demonstrated that this sample contained Banna virus (BAV) belonging to the family Reoviridae in the genus Seadornavirus. BAV has 12 segments of double-stranded viral RNA, and our NGS result detected all 12 segments of BAV, showing 94.3–98.3% nucleotide identity to the previously reported BAV [1] (Supplementary Fig. 2).

We next designed BAV-specific primer set, BANV-F and BANV-R (Table 1), for detection of BAV VP1 gene (Supplementary Fig. 2) and performed RT-PCR. RNA samples from Vero (40 pools), BHK-21 (43 pools), and C6/36 (65 pools with CPEs) cells for BAV were analyzed. Eight samples that showed CPEs in C6/36 cells were positive, but no samples in Vero and BHK-21 cells were positive. BAV was detected from the pools of female Ae. albopictus (4 pools), Ae. aegypti (2 pools), and Cx. tritaeniorhynchus (2 pools) mosquitoes, all of which were captured in Bogor, West Java province (site no. 4 in Fig. 1 and Table 3).

Detection and isolation of permutotetravirus

Next, we selected the pool that showed CPEs in C6/36 cells (pool no. 34 from Ae. albopictus collected in Bogor, 2017) (Supplementary Fig. 1) and was negative by RT-PCR including BAV-specific primers for NGS by MiSeq (Illumina). The sample contained viral RNAs with high homology to Sarawak virus (SWKV) [49], unassigned member of tetravirus. SWKV possesses approximately 5.3 kb of single-stranded positive-sense RNA genome containing two large open reading frames (ORFs) [49]. We detected nearly complete genome sequences (5319 nt, without RACE analysis), showing 88.9% nucleotide identity to that of SWKV strain SWK-M26 (GenBank accession no. NC_040540.1, Supplementary Fig. 3). We tentatively named this isolate Bogor virus (BGV) strain BGV/MQ/5/Bogor/2017.

We designed the BGV-specific primer set, BGV-F and BGV-R (Table 1), and performed RT-PCR for detection of BGV (Supplementary Fig. 3). Because members of the family Permutotetraviridae are known ISVs [10], supernatant from C6/36 cells (65 pools that showed CPEs) were analyzed by RT-PCR. Thirty-four samples were positive for BGV, which were obtained from the pools of Ae. aegypti (5 pools of female), Ae. albopictus (7 pools of female and 9 pools of male), Ar. subalbatus (6 pools of female), Cx. quinquefasciatus (6 pools of female), and Cx. fuscocephala (1 pool of female) mosquitoes collected in West Java (Bogor), Banten, and Central Kalimantan provinces (site nos. 2, 4, and 7 in Fig. 1 and Table 3).

Detection and isolation of alphamesoniviruses

Lastly, we analyzed one pool that showed CPE in C6/36 cells (pool no. 69 from Cx. tritaeniorhynchus collected in Bogor, 2017) (Supplementary Fig. 1) and was negative by RT-PCR for flaviviruses, BAV, and BGV, by using MinION. We found that 4 of 69 reads obtained by MinION showed high homology to Dak Nong virus [27] (one strain of the viral species Alphamesonivirus 3, genus Alphamesonivirus, family Mesoniviridae; GenBank accession no. NC_038297.1) (Supplementary Fig. 4). No viral sequences, except for alphamesoniviruses, were detected in the 69 reads.

We designed alphamesonivirus-specific primers, AMSV-F and AMSV-R (Table 1), and performed RT-PCR for detection of viral genome (Supplementary Fig. 4). Because viruses of the family Mesoniviridae are also known to infect only insect cells [55], RNA samples from C6/36 cells (65 pools with CPEs) were analyzed. Twenty-three samples were positive and originated from Ae. albopictus (2 pools of female), Ar. subalbatus (1 pool of female), Cx. quinquefasciatus (2 pools of female and 1 pool of male), and Cx. tritaeniorhynchus (17 pools of female) mosquitoes in West Java and Banten provinces (site nos. 1, 2, and 4 in Fig. 1 and Table 3). Notably, flaviviruses, BAV, BGV, and alphamesoniviruses were not detected from the same samples.

We successfully isolated 8 insect flaviviruses, 8 BAV, 34 BGV, and 23 alphamesoniviruses from mosquitoes in Indonesia (Table 3). Next, we analyzed phylogenetic relationships of each viral isolate with related viruses.

To date, more than 25 species of insect flaviviruses have been reported globally [3, 6]. We conducted phylogenetic analysis of our eight insect flaviviruses to infer the relationship with the reported insect flaviviruses based on partial NS5 gene sequences (203 nt). As mentioned above, these eight isolates were classified into three species of insect flaviviruses: six CFAV, one AEFV, and one CxFV strain. The phylogenetic tree showed that our isolated flaviviruses formed clusters with each viral species (Fig. 2). Among our six CFAV strains, there were only two synonymous substitutions in the 203 bases, and they formed one subcluster differing from another CFAV from Thailand, Australia, USA, Puerto Rico, and Mexico (Fig. 2). AEFV has been reported from Japan, Italy, USA, Thailand [3], Turkey [46] and China [13]. The result of phylogenetic analysis indicated that our Indonesian AEFV strain AEFV/MQ/29/Bogor/2017 was located as a root of the other AEFV strains (Fig. 2). CxFV was first isolated from Japan [20] and has been reported from many countries, including China, Taiwan, Indonesia, Guatemala, USA, Mexico, Trinidad, Uganda, Argentina, and Brazil, indicating that CxFV is ubiquitous in the world [3]. The phylogenetic analysis indicated that our CxFV strain CxFV/MQ/174/Bogor/2017 was most closely related to CxFV strain Surabaya-2 isolated in 2004 from mosquitoes in Java, Indonesia [22] (Fig. 2).

Fig. 2.

Phylogenetic tree of representative insect-specific flaviviruses based on the partial NS5 gene (203 nt) by the maximum likelihood algorithm. Each branch was assessed by the bootstrap method with 100 replicates under the generalized time reversible model. Sequence data for reference isolates were obtained from the NCBI database, and the country and year of isolation are provided in parentheses. Bold and underlined letters show isolates detected in this study. DDBJ accession numbers are shown in brackets.

Because AEFV was first isolated in Indonesia, we further analyzed complete genome sequences of AEFV/MQ/29/Bogor/2017 by NGS using MiSeq (Illumina) and the analysis were performed using the same method for BGV/MQ/5/Bogor/2017. We successfully determined almost complete genome sequences for AEFV/MQ/29/Bogor/2017 (11,072 nt, without RACE analysis) (GenBank accession no. LC536088). To date, four AEFV strains from Japan (strain Narita-21; NC_012932.1) [21], Thailand (strain Bangkok; KJ741266.1) [5], Turkey (strain KRD32; MK251047.1) [46], and USA (strain SPFLD-MO-2011-MP6; KC181923.1) [19] have been deposited in the NCBI database with full or nearly-full genome sequences (as of December 2019). Nucleotide sequence identity for the polyprotein region (10,026 nt) between AEFV/MQ/29/Bogor/2017 and the other AEFV strains ranged from 89.5% with strain Narita-21 to 88.2% with strain Bangkok. We conducted phylogenetic analysis of the five AEFV strains based on nucleotide sequences encoding the polyprotein region, which indicated that the five AEFV strains divided into three subbranches: Indonesian strain, Thailand strain, and the other strains (Supplementary Fig. 5).

BAV was first isolated from the cerebrospinal fluid of encephalitis patients in China [59]. Subsequently, the virus has been isolated from mosquito specimens collected in Indonesia [7], Vietnam [42], and several provinces of China [34] and also Culicoides midge in China [35, 52]. In this study, we determined the sequence of RT-PCR products for eight BAV isolates and then analyzed phylogenetic relationships with other BAV strains based on the partial VP1 gene (736 nt, Supplementary Fig. 2). BAV were divided into some subclusters, and all our eight BAV strains belonged to the same cluster with Indonesian BAV isolated during 1980–1981 (Fig. 3). There was no nucleotide substitution among the eight BAV isolates obtained in this study.

Fig. 3.

Phylogenetic tree of Banna viruses based on the partial VP1 gene (736 nt) by the maximum likelihood algorithm. Each branch was assessed by the bootstrap method with 100 replicates under the generalized time reversible model. Sequence data for reference isolates were obtained from the NCBI database, and the country and year of isolation are provided in parentheses. Bold and underlined letters show isolates detected in this study. DDBJ accession numbers are shown in brackets.

On the website of the International Committee on Taxonomy of Viruses (ICTV; https://talk.ictvonline.org/), two viral species, Euprosterna elaeasa virus and Thosea asigna virus, are members of the genus Alphapermutotetravirus in the family Permutotetraviridae (assigned one genus for this family) as lepidopteran insect viruses. In addition, several related viruses, such as Newfield virus [57], SWKV [49], Shinobi tetravirus (SHTV) [15], and Egaro virus (no article available, GenBank accession no. NC_030845.1), have been reported from dipteran specimens recently. In this study, we selected 20 BGV isolates from 34 RT-PCR samples positive for BGV and determined their sequences (393 nt, Supplementary Fig. 3). The analysis revealed 23 substitutions causing two amino acid differences. Our isolates showed 87.2–88.5% nt and 95.4–96.1% aa sequence identity to that of SWKV strain SWK-MS26 (NC_040540.1), and 84.2–85.4% nt and 93.1–93.8% aa sequence identity to that of SHTV strain Shinobi (LC270813.1). Phylogenetic analysis for permutotetraviruses based on the 131 aa sequences of the putative RdRp region indicated that our 20 BGV isolates formed one cluster with SWKV and STHV (Fig. 4).

Fig. 4.

Phylogenetic tree of Bogor viruses with the related permutotetraviruses based on putative amino acid sequences (131 aa) from partial open reading frame 1 (ORF1) sequences by the maximum likelihood algorithm. Each branch was assessed by the bootstrap method with 100 replicates under the JTT model. Sequence data for reference isolates were obtained from the NCBI database, and the country and year of isolation are provided in parentheses. Bold and underlined letters show isolates detected in this study. DDBJ accession numbers are shown in brackets.

On the ICTV website, members of the family Mesoniviridae (order Nidovirales, suborder Mesonidovirineae) currently consist of nine species (Alphamesonivirus 1–9). We carried out sequence analysis of 23 RT-PCR products by using the AMSV-specific primer set. We found that 7 of 23 RT-PCR products were overlapped sequences, possibly due to superinfection of 2 or more alphamesoniviruses. Therefore, we used 16 mesonivirus sequences for phylogenetic analysis (336 nt, Fig. 5). The result revealed that they were divided into 2 alphamesonivirus species, 14 isolates of Alphamesonivirus 2 and 2 isolates of Alphamesonivirus 3 (Fig. 5). Alphamesonivirus 2 was obtained from pools of Ae. albopictus, Ar. subalbatus, Cx. quinquefasciatus, and Cx. tritaeniorhynchus mosquitoes, and Alphamesonivirus 3 was obtained from pools of Cx. tritaeniorhynchus mosquitoes (Table 3). There were no substitutions among the 14 isolates of Alphamesonivirus 2 and 2 isolates of Alphamesonivirus 3. We concluded that at least two alphamesonivirus species (Alphamesonivirus 2 and Alphamesonivirus 3), including the possibility for superinfection with another alphamesonivirus species, were distributed in the wild mosquito population in Indonesia.

Fig. 5.

Phylogenetic tree for alphamesoniviruses based on the partial open reading frame 1 (ORF1) gene (312 nt) by the maximum likelihood algorithm. Each branch was assessed by the bootstrap method with 100 replicates under the generalized time reversible model. Sequence data for reference isolates were obtained from the NCBI database, and the country and year of isolation are provided in parentheses. Bold and underlined letters show isolates detected in this study. DDBJ accession numbers are shown in brackets.

DISCUSSION

In this study, we conducted surveillance of mosquito-borne viruses in four islands of Indonesia, Sumatra, Java, Bali, and Borneo, during 2016–2018. A total of 10,015 mosquitoes were collected during the study and used for virus isolation, resulting in isolation of many mosquito-borne viruses belonging to families Flaviviridae (three insect flaviviruses: CFAV, AEFV, and CxFV), Reoviridae (BAV), Permutotetraviridae (tentatively named Bogor virus), and Mesoniviridae (alphamesoniviruses 2 and 3). These findings are valuable for quick identification and diagnosis of viruses in captured mosquitoes in Indonesia and contribute to the understanding of the biodiversity of mosquito-virus relationships and virus taxonomy in tropical regions.

Many viruses in this study were predominantly ISVs except for BAV, and no arboviruses, such as JEV and DENV, were isolated. Previous studies reported that DENV genomes were detected by RT-PCR in 0.12% (36/29,252) [47] and 0.66% (110/16,605) [41] in mosquito population. In a DENV-endemic area in Brazil, the MIR of DENV by RT-PCR in wild Aedes mosquitoes was 19.8 [37]. In this study, we collected approximately 3,000 Aedes mosquitoes (both male and female of Ae. aegypti and Ae. albopictus) (Table 2) and used them for virus isolation, but not for RT-PCR analysis. Therefore, the number of our collected Aedes mosquitoes may be insufficient for detection of Aedes-borne arboviruses by virus isolation. Survey of JEV in Culex mosquitoes revealed that the MIR of JEV in Indonesia was 0.06 (1/1,485) [16]. Recently in China, 10 of 511 pools of Culex mosquitoes were positive by RT-PCR for detection of JEV genome, and the MIR of JEV was estimated at 0.91 [14]. Thenmozhi et al. reported that JEV-antigen was detected in 4 of 580 mosquitoes in India by enzyme-linked immunosorbent assay (ELISA) and the MIR was estimated at 6.9 [53]. Because of the dominant religion in Indonesia, rearing of domestic pigs has been restricted except for in some provinces such as Bali [23]. In this study, we collected 1,030 adult female Cx. tritaeniorhynchus mosquitoes, known as the main vector for JEV, mainly in Bogor, West Java province (Table 2), where domestic pigs are rare [39]. Therefore, the number of collected Cx. tritaeniorhynchus mosquitoes may not be adequate for isolation of JEV from mosquitoes. We should improve the method for efficient collection of arbovirus-infected mosquitoes for further entomological surveillance and risk assessment of arboviral diseases.

Three insect flaviviruses, CFAV, AEFV, and CxFV, were isolated from Ae. aegypti, Ae. albopictus, and Cx. quinquefasciatus, respectively (Fig. 2 and Table 3). Although CFAV and CxFV have been reported from Surabaya city, Java, in Indonesia [20, 21], the isolation of AEFV is the first report in Indonesia. In the phylogenetic analysis, our CxFV isolate CxFV/MQ/174/Bogor/2017 formed a subcluster with CxFV strain Surabaya-2 (Fig. 2). Partial genome sequences encoding NS5 gene of Indonesian CFAV strain Surabaya-10 [21] was deposited in GenBank (accession no. AB488428.1) and showed overlap of only 91 bp with sequences of our CFAV strain. Additionally, the highest homology between Surabaya-10 and our CFAV strain was observed in this 91-bp region (data not shown). These results implied that insect flaviviruses have been maintained and circulate in Java, Indonesia.

BAVs were isolated from several species of female mosquitoes, Ae. aegypti, Ae. albopictus, and Cx. tritaeniorhynchus (Table 3). A previous study reported that BAV was also isolated from Cx. pseudovishnui, Ae. vagus, and Cx. pipiens at several sites in Indonesia during 1980–1981 [7]. Other mosquito species, such as Anopheles sinensis and Ae. vexans, and biting midge (Culicoides sp.) have been also reported as potential vectors for BAV in China [34, 52, 58]. BAV was first isolated from cerebral fluid of human patients with encephalitis [59], and several seroepidemiological studies on BAV infection in human and animals have suggested the potential of BAV to cause arboviral diseases [52]. These findings suggested that BAV has been maintained among several mosquito species in various Asian countries and probably circulate between mosquitoes and vertebrate hosts in nature [35]. However, we did not detect any BAV isolates from mosquito pools inoculated on both mammalian Vero and BHK-21 cell lines. Further investigation of transmission dynamics and pathogenicity of BAV strains in vertebrates is required.

On ICTV taxonomy, the viral family Permutotetraviridae currently consists of one genus Alphapermutotetravirus and two species, Euprosterna elaeasa virus and Thosea asigna virus, are assigned. They are known as lepidopteran insect viruses. In 2010, the viral genome of permutotetra virus-like virus, named Newfield virus, was discovered from wild-collected Drosophila flies [57]. In 2017, SWKV was isolated from a pool of male Ae. albopictus mosquitoes in Sarawak State (northwest Borneo) of Malaysia as a new member of tetravirus [49]. Our colleagues reported the isolation of SHTV strain Shinobi from the Ae. albopictus mosquito cell line, JCRB, in 2018 as a new member of the family Permutotetraviridae [15]. Overall, three new permutotetraviruses have been reported as a new group of probable dipteran insect tetraviruses. In this study, we successfully determined the almost full genome sequences of BGV strain BGV/MQ/34/Bogor/2017 (5319 nt without terminal regions), which contains two putative ORFs (Supplementary Fig. 3, GenBank accession no. LC536112). The putative amino acid sequences of BGV ORFs showed high homologies to those of SWKV strain SWK-M26 (GenBank accession no. NC_040540.1; 96.5% for ORF1 and 91.3% for ORF2) and SHTV strain Shinobi (GenBank accession no. YP_009553485.1; 94.8% for ORF1 and 88.5% for ORF2).

We obtained 34 BGV isolates from various mosquito species, Ae. aegypti, Ae. albopictus, Ar. subalbatus, Cx. quinquefasciatus, and Cx. fuscocephala, at several sites in Indonesia (Java and Borneo), suggesting that this virus may be widely distributed in tropical areas and mosquito populations. However, it is unknown whether BGV can transmit among different mosquito species.

Twenty-three alphamesoniviruses were isolated from various mosquito species (Table 3), and phylogenetic analysis of these isolates revealed that they were classified as at least two viral species, Alphamesonivirus 2 and Alphamesonivirus 3 (Fig. 5). Alphamesonivirus 2, alias Karang Sari virus, was isolated from specimens of a pool of Cx. vishnui mosquitoes in 1981 at Karang Sari, Banten province, Indonesia [55]. We isolated Alphamesonivirus 2 from mosquito pools of Ae. albopictus, Ar. subalbatus, Cx. quinquefasciatus, Cx. tritaeniorhynchus, in Java in 2017. By contrast, Alphamesonivirus 3, alias Dak Nong virus, was discovered from a pool of Cx. tritaeniorhynchus mosquitoes in Vietnam [27]. Another related strain of Alphamesonivirus 3, Kamphaeng Phet virus, was found in Thailand, but no information on mosquitoes was available [55]. In this study, we obtained Alphamesonivirus 3 only from pools of Cx. tritaeniorhynchus mosquitoes in Java, suggesting that alphamesoniviruses might possess host-specificity to mosquito species. We found the mixed growth of two or more alphamesonivirus species in the same C6/36 cultures (7 samples), although mixed growth has not been observed among flaviviruses, BAV, and BGV. We hypothesize that no competition occurred among alphamesoniviruses or genetically related viruses in C6/36 cells. These interactions among viruses in mosquitoes should be analyzed to clarify virus-vector interaction.

Mosquito-borne viruses are considered a significant cause of major health problems in Indonesia. Although we did not isolate human-related arboviruses in this study, we believed that this large surveillance of mosquitoes and virus isolation provides a more accurate view of the prevalence of arboviruses. Furthermore, the presence and distribution of ISVs in mosquito populations raise interesting issues about their possibility to interact with pathogenic viruses and to be the next emerging infectious pathogen. Monitoring and long-term surveillance of mosquito-borne viruses are critical to prevent and control emerging and re-emerging arboviral diseases in Indonesia.