De novo profiling of RNA viruses in Anopheles malaria vector mosquitoes from forest ecological zones in Senegal and Cambodia

Abstract

Background: Mosquitoes are colonized by a large but mostly uncharacterized natural virome of RNA viruses, and the composition and distribution of the natural RNA virome may influence the biology and immunity of Anopheles malaria vector populations.

Results: Anopheles mosquitoes were sampled in malaria endemic forest village sites in Senegal and Cambodia, including Anopheles funestus, Anopheles gambiae group sp., and Anopheles coustani in Senegal, and Anopheles hyrcanus group sp., Anopheles maculatus group sp., and Anopheles dirus in Cambodia. The most frequent mosquito species sampled at both study sites are human malaria vectors. Small and long RNA sequences were depleted of mosquito host sequences, de novo assembled and clustered to yield non-redundant contigs longer than 500 nucleotides. Analysis of the assemblies by sequence similarity to known virus families yielded 115 novel virus sequences, and evidence supports a functional status for at least 86 of the novel viral contigs. Important monophyletic virus clades in the Bunyavirales and Mononegavirales orders were found in these Anopheles from Africa and Asia. The remaining non-host RNA assemblies that were unclassified by sequence similarity to known viruses were clustered by small RNA profiles, and 39 high-quality independent contigs strongly matched a pattern of classic RNAi processing of viral replication intermediates, suggesting they are entirely undescribed viruses. One thousand five hundred sixty-six additional high-quality unclassified contigs matched a pattern consistent with Piwi-interacting RNAs (piRNAs), suggesting that strand-biased piRNAs are generated from the natural virome in Anopheles. To functionally query piRNA effect, we analyzed piRNA expression in Anopheles coluzzii after infection with O'nyong nyong virus (family Togaviridae), and identified two piRNAs that appear to display specifically altered abundance upon arbovirus infection.

Conclusions: Anopheles vectors of human malaria in Africa and Asia are ubiquitously colonized by RNA viruses, some of which are monophyletic but clearly diverged from other arthropod viruses. The interplay between small RNA pathways, immunity, and the virome may represent part of the homeostatic mechanism maintaining virome members in a commensal or nonpathogenic state, and could potentially influence vector competence.

Keywords: Anopheles; Insect specific virus; Malaria vector; RNA virus; Virome; Virus genome assembly.

Fig. 1

Taxonomic profile of Anopheles sample pools. Relative abundance values of Anopheles species were computed by mapping long RNAseq reads to mitochondrial cytochrome C oxidase subunit I gene sequences from the Barcode of Life COI-5P Database. Taxa represented by > 100 sequence reads and 1% frequency in the sample pool were plotted in pie charts. White wedges in pie charts represent the combined proportion of all sequence matches that were individually present at less than 1% frequency in the sample. All data are presented in tabular form in Additional file 1: Table S1

Fig. 2

Phylogenetic tree of reference and novel virus assemblies from the Bunyavirales order. Maximum-likelihood phylogeny based on RNA-dependent RNA polymerase (RdRP) predicted peptide sequences of viruses from the Bunyavirales order. Novel viruses characterized in the current study (red name labels) are placed with reference viruses (black name labels) within the Phasmavirus clade and in a basal position of the Phlebovirus-Tenuivirus clade. Node robustness is indicated by bootstrap values (number of replicates supporting the node), indicated by color of the dot at the branch base, see key. Protein lengths and functional status of RdRP peptide sequences from novel viruses in the current study are included to distinguish between complete and partial and/or non-functional pseudogenes (indicated by label “pseudogenized”, functional status also shown in Additional file 2: Table S2 and Additional file 3: Table S3). Average protein size of reference virus RdRP genes is 2496 amino acids

Fig. 3

Phylogenetic tree of reference and novel virus assemblies from the Mononegavirales order. a Maximum-likelihood phylogeny based on RNA-dependent RNA polymerase (RdRP) predicted peptide sequences of viruses from Mononegavirales order. Novel virus assemblies characterized from Cambodia and Senegal Anopheles samples (red name labels) are placed with reference viruses (black name labels), predominantly within the Dimarhabdovirus clade and as close relative of the Nyamivirus clade. Node robustness is indicated by bootstrap values (number of replicates supporting the node), indicated by color of the dot at the branch base, see key. Protein lengths and functional status of RdRP peptide sequences from novel viruses in the current study are included to distinguish between complete and partial and/or non-functional pseudogenes (indicated by label “pseudo”, functional status indicated in Additional file 2: Table S2 and Additional file 3: Table S3). Average protein size of reference virus RdRP genes is 2098 amino acids. b Genome comparison of novel and reference Xincheng Mosquito Viruses, which are too diverged to align at the nucleic acid sequence level. Grey blocks represent peptide sequence homology regions between compared sequences. The nucleotide sequences of the entire viral contigs, and not only of the RdRP gene as in (a), were translated and used to search the translated nucleotide database with TBLASTX. The viruses display recognizable relatedness over their genomes, despite geographic distance and nucleotide sequence divergence. Color intensity indicates identity levels from TBLASTX results (values indicated in key)

Fig. 4

Virus abundance profiles across mosquito sample pools based on long and small RNA sequence mapping. Heatmap of log2-transformed reads per kilobase per million reads (RPKM) abundance values of novel virus assemblies identified from Cambodia and Senegal sample pools based on long and small RNA sequence libraries. Broadly similar viral abundance profiles are detected in sample pools by the long and small RNA sequence data. Representation of particular viruses is uneven among mosquito sample pools, suggesting inter-individual mosquito differences for virus carriage. X-axis, Anopheles sample pools from Cambodia, Cam, and Senegal, Dak; y-axis, names of 115 assembled virus contigs displaying sequence similarity to known virus families (Additional file 2: Table S2 and Additional file 3: Table S3)

Fig. 5

Small RNA size profiles of novel virus assemblies from Cambodia and Senegal sample pools. Hierarchical clustering of 88 novel virus assemblies based on Pearson correlation of small RNA size profiles. The 88 viruses were the members of the 115 novel virus set meeting the threshold of at least 100 small RNA reads mapped to the viral contig, to assure reliable small RNA size profiling. Small RNA reads that mapped to each of the 88 virus assemblies were extracted, and their size distributions were normalized with a z-score transformation. Heatmaps indicate the frequency of small RNA reads of size 15 to 35 nucleotides that map over the positive strand (left panel) and negative strand (right panel) of the reference sequence indicated on the y-axis. The x-axis indicates the size in nucleotides of the small RNAs mapped. Four main clusters were defined (indicated by numbers on the left of each panel) based on these small RNA size profiles. The profile in Cluster 3 is enriched for 21 nucleotide reads mapping over both positive and negative strands, characteristic of the classical small interacting RNA (siRNA) product size profile

Fig. 6

O’nyong nyong arbovirus infection influences expression of candidate piRNA genes in Anopheles coluzzii. Anopheles coluzzii mosquitoes were challenged with O’nyong nyong virus (ONNV) by feeding an infectious bloodmeal or an uninfected control bloodmeal, and small RNAs expressed during the primary infection at 3 d post-bloodmeal were sequenced. Analysis using Cuffdiff highlighted two candidate piRNA genes that displayed decreased abundance of mapped small RNAs in ONNV infected samples (see Results, piRNA loci XLOC_012931 and XLOC_012762). Here, the small RNA sequence reads mapping to the two candidate piRNA loci were quantified using the Integrative Genomics Viewer normalized to the library size, and the difference between ONNV infected and uninfected samples tested statistically. X-axis indicates candidate piRNA locus, y-axis indicates percentage of normalized small RNA reads mapping to the piRNA gene. ONNV-infected mosquitoes, red bar; uninfected control mosquitoes, black bar. Experiments were done in two biological replicates, error bars indicate standard deviation. Locus XLOC_012931, Chi-squared = 77.36, df = 1, p-value< 2.2e-16 (ONNV-infected mean mapped reads = 36 ± 141,421,356, mean total reads = 19,193,551 ± 8,555,908.61, ONNV-uninfected mean mapped reads = 160 ± 14,1,421,356, mean total reads = 19,167,336 ± 3,962,902.88052); and locus XLOC_012762, Chi-squared = 75.78, df = 1, p-value< 2.2e-16 (ONNV-infected mean mapped reads = 51 ± 19,09, mean total reads = 19,193,551 ± 8,555,908.61, ONNV-uninfected, mean mapped reads = 184 ± 848,528,137, mean total reads = 19,167,336 ± 3,962,902.88)