Metagenome analysis of viruses associated with Anopheles mosquitoes from Ramu Upazila, Cox’s Bazar District, Bangladesh

Study area

The mosquito collection was carried out in Ramu Sub-district of Cox’s Bazar District, Bangladesh from November 2018 to April 2019. The site is situated between 21°17′ and 21°36′ north latitudes and in between 92°00′ and 92°15′ east longitudes (Fig. 1); it covers approximately 391.91 km² with a population of 266,640 people at a density of 680.7/km² according to a 2011 census. Ramu is bounded by Chakaria and Cox’s Bazar sub-districts on the north, Naikhongchhari and Ukhia sub-districts on the south, Naikhongchhari sub-district on the east, Cox’s Bazar Town and the Bay of Bengal on the west. The landscape comprises paddy fields on the plains and unused shrub lands or teak plantations in hilly areas. Entomological investigations were conducted in the following three villages: Horitola of Rashidnagor union (21°29′817″ N, 92°6′082″ E), VIP Pahar and Forest Office Purbopara of Jowarianala union (21°29′028″ N, 92°6′825″ E) and Moheshkum of Kawerkop union (21°26′062″ N, 92°9′294″ E).

Mosquito collections

Adult active mosquitoes were collected only from human dwellings using battery operated CDC miniature light traps Model 1012 (John W. Hock Company, Gainesville, FL, USA) between hours 1800 to 0600 following the procedure described by Alam et al. (2010). Traps were set both inside and outside the main living room, i.e., indoor and outdoor, respectively. In each month, 18 traps were operated. Six traps (three indoor and three outdoor) were installed per night in each village. Thus, in a 6-month collection period a total of 108 traps were installed at the study sites.

Resting adult mosquitoes were collected by pyrethrum spray catches (PSCs) between hours 1900 to 1000 (Ndiath et al., 2011). A single bedroom in each house was selected and all windows, doors, and other escape passages of that room were sealed. Pyrethrum insecticides (the aerosols generally used for domestic purposes) were sprayed in the room for 9–12 s (by dividing the room into 3–4 parts based on the size of the room) and the door was closed for approximately 15 min. Immobilized mosquitoes were collected using mouth aspirator with HEPA filter and kept in plastic vials over the course of 15–25 min. Windows and doors were re-opened for proper ventilation and all members of the house were requested to remain outside for 15 min. Each month, five PSCs were carried out per village. In total, 15 PSCs were conducted spanning three villages per month over 6 months, for a total of 90 PSCs.

Collected Anopheles mosquitoes were morphologically identified using standard taxonomic keys (Rattanarithikul et al., 2006) and kept in individual vials containing silica gel desiccant for storage at room temperature at the International Centre for Diarrhoeal Disease Research, Bangladesh (icddr,b). The room temperature stored samples were shipped to Walter Reed Army Institute of Research on July 25, 2019, and stored under −80 °C until extraction for nucleic acids purification.

Nucleic acids extraction and purification from mosquitoes

Male Anopheles mosquitoes or female Anopheles mosquitoes with no visible engorged blood meal (unfed female) were sorted into pools of up to ten individuals of same sex and species from the same collection date and household and transferred into a 2-ml screw-cap microcentrifuge tube with one 3-mm glass bead and two 1-mm glass beads (Sigma-Aldrich, Burlington, MA, USA). For each pool, 0.6 ml of cell growth medium, Eagle’s Minimum Essential Medium (Quality Biological, Inc., Gaithersburg, MD, USA) supplemented with 10% fetal bovine serum (Thermo Fisher Scientific, Waltham, MA, USA) was added, followed by homogenization on a Mini-Beadbeater 16 (BioSpec Products, Inc., Bartlesville, OK, USA) for 1 min and then centrifugated at 4 °C at 8,000 × g for 20 min. Prior to extraction, free nucleic acids were first digested by adding 12 µl of enzyme mixture, containing 10.4 units of DNase I (New England BioLabs, Ipswich, MA, USA), 53.6 units of Benzonase (Sigma-Aldrich, Burlington, MA,USA) and 2.24 units of RNase A (Qiagen, Venlo, USA) into each well of a 96-well KingFisher deep-well plate (Thermo Fisher Scientific, Waltham, MA, USA). For each sample, 188 µL of clear supernatant was transferred into the deep-well plate, incubated at 37 °C for 90 min and then proceeded to nucleic acid extraction on KingFisher Flex magnetic bead purification system using MagMAX Viral/Pathogen (MVP) DNA/NA kit (Thermo Fisher Scientific, Waltham, MA, USA) per manufacturer’s protocol.

Random amplification, metagenome sequencing and data analysis

Unbiased reverse transcription and PCR cDNA amplification was conducted as described previously (Sanborn et al., 2021; Hang et al., 2012). In brief 9.7 µl of purified RNA was pre-treated with DNase I and then annealed with random hexamer primers at 65 °C for 5 min. SuperScript III Reverse Transcriptase (Thermo Fisher Scientific, Waltham, MA, USA) and Platinum Taq DNA Polymerase High Fidelity (Thermo Fisher Scientific, Waltham, MA, USA) were used in a reverse transcription and PCR amplification, respectively. The thermal cycler program for PCR was: 94 °C for 3 min to activate the Taq DNA polymerase followed by 35 cycles of amplification at 94 °C for 30 s, 55 °C for 30 s and 72 °C for 1 min. The random amplicons were purified using AMPure XP beads (Beckman Coulter, Brea, CA, USA) and quantified using Quant-iT PicoGreen dsDNA Assay Kit (Thermo Fisher Scientific, Waltham, MA, USA). A NGS library was prepared using 25 ng of random amplicons and the Illumina DNA Sample Prep kit (Illumina, San Diego, CA, USA) following the manufacturer’s protocol. The indexed libraries were quantitated using PicoGreen dsDNA assay and pooled at equal concentrations. Concentration and size distribution for the final library pool was measured using a Tapestation instrument and Agilent DNA 5,000 assay kit (Agilent Technologies, Santa Clara, California, USA). The libraries were sequenced using Illumina MiSeq NGS system and Miseq Reagent Kit v3 (600 cycles) in a 2 × 300 bp paired-end run.

Illumina MiSeq sequence read data were analyzed using an in-house pathogen discovery pipeline (Kilianski et al., 2015), which consists of quality processing, sequence assembly, nucleotide BLAST analyses of both assembled contigs and all remaining unassembled reads against the NCBI GenBank non-redundant nucleotide (nr/nt) database, and taxonomic identification of organisms in each sample. The Chan Zuckerberg ID Short-Read MNGS workflow (https://czid.org/; https://github.com/chanzuckerberg/idseq-workflows), a cloud-based metagenomics platform (Kalantar et al., 2020), was also used for visualization of mNGS results of selected samples. The sequence contigs with a length of 1 kb or greater were retained for further analysis. MAFFT version v7.475 with the L-INS-I strategy was used to align nucleotide contigs and publicly available background sequences (Katoh & Toh, 2008). IQ-TREE version 1.6.12 (Nguyen et al., 2015) with 1,000 ultrafast boostrap replicates (“--bb 1,000”) was used in phylogenetic analysis. IQ-TREE ModelFinder was used to determine the best-fit model. FigTree (http://tree.bio.ed.ac.uk/software/figtree/) v1.4.4 was used for phylogenetic tree visualization and figure generation.

Endogenous viral element (EVE) classification

To assess the shared ancestry of, and screen for potential EVEs within, the 64-contig dataset with sequence lengths >1 kbp, we first retrieved the aligned nucleotide window of all associated BLASTn hits to any viral target with an alignment length ≥ 200 bp and a bitscore ≥100 from the NCBI nt database. Additionally, we conducted a second BLASTn search using our 64 queries against a database comprising all Anopheles mosquito genomes available at the time of analysis (n = 89; as of 21 November 2023) available in the NCBI whole-genome shotgun (WGS) database and extracted any hits as above. Each set of extracted hits (viral and mosquito) was clustered at 99% nucleotide identity using cd-hit-est (Fu et al., 2012) to remove database sequences that were returned in multiplicity and thus redundant due to closely related query sequences and overlapping BLAST homology windows. These viral and mosquito-derived representative sequences were then aligned together with the 64 viral contig queries using mafft-einsi (Katoh & Toh, 2008) and a maximum-likelihood phylogeny was constructed from the alignment using IQTREE2 under automatic model selection with nodal supports assessed using 2,000 rapid bootstrap approximations (Nguyen et al., 2015).

Metagenomic sequencing of Bangladesh Anopheles virome

In total, 150 pools (644 individual mosquitoes) of Anopheles mosquitoes were sequenced; An. vagus comprised over half of the total specimens (54%), with an additional 14 other Anopheles species (42%) and a small number of unidentified An. spp. (4.1%) available for mNGS analysis (Table S1). An. gambiae and An. stephensi mosquitoes were used in concomitant malaria studies and not analyzed in this study. Males comprised 63 pools (293 individuals) and females, 87 pools (356 individuals) of un-engorged mosquitoes. A total of 123.4 million of MiSeq sequence reads were obtained for metagenomic analysis. Among the de novo assembled contigs with 500 or more mapped sequence reads, 107 contigs were identified as viral via BLASTn analysis against the NCBI nr/nt database. The identified virus sequences were found in 53 of the 150 pools (35.3%). The identified viruses were from 12 established viral families, including Dicistroviridae, Flaviviridae, Hantaviridae, Iridoviridae, Narnaviridae, Partitiviridae, Peribunyaviridae, Phasmaviridae, Rhabdoviridae, Sedoreoviridae, Solemoviridae, Totiviridae, and unclassified Bunyavirales, unclassified Riboviria, or unclassified viruses (Fauver et al., 2016; Piegu et al., 2014; Coffey et al., 2014; Truong Nguyen et al., 2022). These 64 viral scaffolds from 40 Anopheles pools with lengths of 1 kb or greater and manually curated (Table 1) were retained for further analysis.

Anopheles associated viruses in Ramu, Bangladesh

The viral contigs summed to 147,106 bp, with lengths ranging from 6,957 bp to 1,012 bp, with an average length of 2,263 bp and mean size of 1,590 bp. The BLASTn alignment showed that many sequences exhibited low nucleotide sequence identity (under 75%) with known sequences in aligned regions (Table 1). Annotations derived from BLASTn homology indicate these 68 protein coding sequences (CDS) comprise 6 nucleocapsid, 16 hypothetical protein, 19 glycoprotein, and 27 RNA-dependent RNA polymerase genes. The virus sequences obtained in this study are named after the sample collection location, the Ramu subdistrict, and based on sequence similarities with known virus sequences in GenBank (Table 1). The partial sequence for Ramu Anopheles Seadornavirus strain KPSC3P0946 (family Sedoreoviridae, genus Seadornavirus), 2,435 bp in length encoding RNA-dependent RNA polymerase (RdRp), was most related with segment 1 of Mangshi virus strain DH13M041 (KR349187) (Wang et al., 2015), isolated from Culex tritaeniorhynchus in Yunan Province, China, with nucleotide identity of 85.7%. The sequence shares 66.1% nucleotide identity with Banna virus (NC_004211.1) which can infect human, animal and insects (Attoui et al., 2000). Ramu Anopheles Sobemo-like virus strain JLT1P0124, partial genome sequence of 2,678 bp, showed low nucleotide sequence identity of 65% with viruses from family Solemoviridae, which were found associated with arthropods, plants, and possibly bats (Somera et al., 2021).

Metagenomic analysis identified sequences of Orthophasmavirus (order Bunyavirales, family Phasmaviridae) in 24 (65%) out of the 40 pools and sequences of rhabdovirus (order Mononegavirales, family Rhabdoviridae) in 29 pools (72.5%) (Fig. 2). Specifically, the reference genome sequences for RNA-dependent RNA polymerase (L segment, NC_031307), glycoprotein (M segment, NC_031308), and nucleocapsid (S segments, NC_031310) of Wuhan mosquito orthophasmavirus 1 (genus Orthophasmavirus, species Wuhan mosquito orthophasmavirus 1 or Orthophasmavirus wuhanense) shared homology with viruses identified in this study, with nucleotide identities of approximately 70% for L segment and 65% for M segments. Alignment of M segment nucleotide sequences of Ramu Anopheles Orthophasmaviruses from this study revealed a cluster of 14 closely related strains (Fig. 3), with nucleotide identity over 97% among the strains.

Thirteen rhabdovirus sequences were identified and designated as Ramu Anopheles Rhabdovirus (Fig. 4). All genome sequences, except for strain KLT5P0681, are nearly identical to recently reported Cambodia Anopheles Rhabdovirus (family Rhabdoviridae) (OR479699–OR479701) (Mohamed Ali et al., 2023), with nucleotide identities of >98% and suggests a high prevalence of related Anopheles rhabdovirus in Southeast Asia. The Ramu Anopheles Rhabdovirus strain KLT5P0681, partial genome sequence of 1,548 nt, has low homology with all other rhabdoviruses from this study and GenBank, with nucleotide identity of 65–70% in aligned regions and bitscores under 200. It shares nucleotide identity of 74.1% and alignment score of 413 with sequence LP1.2_CulexRhabdovCamb, which was identified in Anopheles mosquitoes in Cambodia (Mohamed Ali et al., 2023). Other previously recognized rhabdovirus sequences that were identified in our phylogenetic analysis include Beaumont viruses from Anopheles mosquitoes in Cambodia (Mohamed Ali et al., 2023), Beaumont virus strain 6 from Anopheles mosquitoes in Australia (Coffey et al., 2014), Xanthi rhabdovirus isolate RC1 and Evros rhabdovirus 2 isolate RB2 from Anopheles mosquitoes in Greece, and Culex pseudovishnui rhabdo-like viruses from Japan (Faizah et al., 2020).

Hubei virga-like viruses (clade Riboviria) are a group of unclassified viruses, which contain protein conserved domain cd23251, an RNA-dependent RNA polymerase (RdRp) in the family Virgaviridae of positive-sense single-stranded RNA [(+)ssRNA] viruses. Sequences that are closely related to Hubei virga-like virus sequences were identified in five Anopheles pools (Fig. 5). Four Ramu Anopheles Virgalike viruses share significant sequence similarity with Hubei virga-like virus 21 strain mosHB236537 (NC_033192, associated with unspecified mosquito species in China) (Shi et al., 2016) and isolate XY8801 (OL700070, Anopheles in China), with amino acid sequence identity of about 80–85%. Sequence for Ramu Anopheles Virgalike virus strain JLT3P1702 is more closely related to Hubei virga-like virus 23 isolate XY43688 (OL700071, Anopheles in China), which is phylogenetically distinct from strain 21. Besides, neither Hubei virga-like viruses 21 and 23 nor Ramu Anopheles Virgalike virus sequence from this study has significant nucleotide identity with any other virga-like viruses identified in Culex mosquitoes or other insects.

Virus sequences identified in male Anopheles mosquitoes

Ten of the 50 virus sequences (20%), which have assembled contigs greater than 1 Kb, were present in male Anopheles (Table 1). These include five orthophasmaviruses, three rhabdoviruses, one sobemolike virus, and one virgalike virus (Table 1), all of which were also seen in female Anopheles in this study and other reports (Figs. 3–5).

Endogenous viral elements

Identification of endogenous viral elements (EVEs) is crucial for future metaviromic work, to avoid misidentification as actively replicating viruses. We identified at least 47 out of the 64 filtered Anopheles viral contigs in clades that shared at least one homologous EVE identified in an Anopheles genome assembly using BLAST searches (Table 2) coupled with phylogenetic analysis (Fig. S1) to confirm monophyly. The largest group of such elements contained 17 of our de novo contigs (12 derived from female and five from male Anopheles mosquitoes), each with homology to the Orthophasmavirus strain Wuhan mosquito virus 1 glycoprotein and integrated as EVEs in the genomes of Anopheles epiroticus, An. cracens, and An. maculatus (Fig. S1 clade A). Remaining homologous elements were identified as derived from: (1) virga-like virus (one contig with top BLAST hits to an element encoded on the genomes of An. atroparvus and An. farauti and six with homology to an RdRP-encoding element currently reported only within An. farauti (Fig. S1 clades B & C, respectively); (2) two nucleocapsid-encoding elements with homology to Anopheles Orthophasmavirus shared broadly with An. maculatus, An. stephensi, An. epiroticus, An. minimus, An. moucheti, An. farauti, An. bellator, An. darlingi, An. ziemanni, An. merus and An. cracens (Fig. S1 clade D); the integration of this EVE may thus pre-date the divergence of major Anopheles species; (3) three contigs encoding nucleocapsid proteins with homology to Anopheles Bunyavirus and present in the genomes of An. koliensis, An. sinensis, An. ziemanni, An. coustani and An. nili (Fig. S1, clade E).; (4) a clade containing twelve RdRP-encoding contigs with homology to Cambodia Anopheles Rhabdovirus with evidence of EVE integration in the genomes of An. rivulorum, An. epiroticus, An. maculatus, and An. funestus (Fig. S1 clade F) that derives from a lineage also containing a separate and distinct Rhabdovirus and associated EVE present in An. merus, An. bwambae, An. quadriannulatus and An. moucheti (Fig. S1 clade G); (5) a single contig encoding an RdRP with homology to Anopheles totivirus and broadly present in the genomes of An. cracens, An. farauti, An. maculatus, An. epiroticus, An. marshallii, An. nili, An. vaneedeni, An. funestus, An. paraensis, An. sinensis, An. atroparvus, and An. punctulatus (Fig. S1 clade H). The number of EVEs derived from viruses related to those described here is evidence for ancient and long-standing evolutionary relationships between these viruses and Anopheles mosquitoes.